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Abstract 

The monotonicity and stability of difference schemes for, in general, 
hyperbolic systems of conservation laws with source terms are studied. 
The basic approach is to investigate the stability and monotonicity of a 
non-linear scheme in terms of its corresponding scheme in variations. Such 
an approach leads to application of the stability theory for linear equation 
systems to establish stability of the corresponding non-linear scheme. The 
main methodological innovation is the theorems establishing the notion 
that a non-linear scheme is stable (and monotone) if the corresponding 
scheme in variations is stable (and, respectively, monotone). Criteria are 
developed for monotonicity and stability of difference schemes associated 
with the numerical analysis of systems of partial differential equations. 
The theorem of Friedrichs (1954) is generalized to be applicable to varia- 
tional schemes with non-symmetric matrices. A new modification of the 
central Lax-Friedrichs (LxF) scheme is developed to be of the second order 
accuracy. A monotone piecewise cubic interpolation is used in the cen- 
tral schemes to give an accurate approximation for the model in question. 
The stability and monotonicity of the modified scheme are investigated. 
Some versions of the modified scheme are tested on several conservation 
laws, and the scheme is found to be accurate and robust. As applied 
to hyperbolic conservation laws with, in general, stiff source terms, it is 
constructed a second order scheme based on operator-splitting techniques. 



1 Introduction 



We are mainly concerned with the stabihty and monotonicity [HJ of difference 
schemes for hyperbohc systems of conservation laws with source terms. Such 
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systems are used to describe many physical problems of great practical impor- 
tance in magneto-hydrodynamics, kinetic theory of rarefied gases, linear and 
nonlinear waves, viscoelasticity, multi-phase flows and phase transitions, shal- 
low waters, etc. (see, e.g., [7], [H], [S], [3D], [35], [SZj, [S], 02], US, [S], 
We will consider a 1-D system of the conservation laws written in the following 
form 

l^ + ^f (U) = iq(u), X e M, < i < Tn,ax; u(x,Olt=o = U°(2^)' (1) 

where u — {ui, U2, . . . , um}'^ is a vector-valued function from R x [0, +c>o) into 
an open subset C M^^, f (u) = {/i (u) , /2 (u) , . . . , /m (u)}"^ is a smooth 
function (flux- function) from fl^ into M^, q (u) = {qi (u) , (72 (u) , . . . , qm {u)}'^ 
denotes the source term, r > denotes the stiffness parameter, u° (x) is either 
of compact support or periodic. We will assume that r = const without loss 
of generality. We will assume that the system ([1]) is hyperbolic and, hence, the 
Jacobian matrix of f (u) possesses M linearly independent eigenvectors (see, 
e.g., [U]). In addition, it is assumed in this paper that 

sup IIAII2 < Amax < 00, A = , (2) 

uen^ o^u 

and all eigenvalues, ^j, ~ ^j, (u) , of the Jacobian matrix G (=9q (u) /9u) have 
non-positive real parts, i.e. 

Re^k (u) < 0, Vfc, Vu ef^u. (3) 
Here and in what follows ||M||p denotes the matrix norm of a matrix M 

induced by the vector norm ||v||p = (^. jwil*')^^'', and ||M|| denotes the matrix 
norm induced by a prescribed vector norm. R and C denote the fields of real 
and complex numbers, respectively, and K denotes either of these fields. 

In the numerical solution of the, in general, stiff (r <C 1) system ([1]), one is 
seeking to establish a numerical scheme that would be robust enough to eradi- 
cate spurious oscillations, i.e. a monotone scheme jlOJ. At the present time there 
are several notions of scheme monotonicity. The notion of 'monotonicity pre- 
serving scheme' originally appeared in Godunov [22| . Such a scheme transforms 
a monotone increasing (or decreasing) function v (x) of space coordinate x at 
a time level t into a monotone increasing (or decreasing, respectively) function 
V (x) at the next time level t + At. Thus, with the use of a monotonicity preserv- 
ing scheme, any discontinuity in the initial monotone data can be smeared in 
succeeding time steps but cannot become oscillatory. Nowadays monotonicity 
preserving schemes are also known as, e.g., monotonicity conserving iterations 
(or methods) [29], monotone schemes (e.g., [1], [34], [44]), monotonicity pre- 
serving methods [40], and Godunov-monotone schemes [10]. Harten et al. [25] 
put forward their own definition of scheme monotonicity as follows: a differ- 
ence scheme, Vi = H {vi-k, Vi-k+i, ■ ■ ■ , Vi+k), is said to be monotone if is a 
monotone increasing function of each of its arguments. The following defini- 
tion is due to Samarskiy: a scheme is regarded as monotone if the boundary 
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maximum principle is maintained [3S] (see also, e.g., [321 P- 183], [TU], [H])- A 

difference scheme may also be referred to as monotone if a maximum princi- 
ple, e.g., the boundary maximum principle, the region maximum principle, the 
maximum principle for inverse column entries, the maximum principle for the 
absolute values, etc', holds for this scheme [TT]. A further important notion 
of difference scheme monotonicity was, in fact, done in [5D] (see also [T^). A 
scheme will be referred to as monotone if it is monotonicity preserving [22] and 
transforms a "A-function" (or "V-function" ) into a "/z - function" (or into an 
"77 - function", respectively). Here and in what follows, a scalar grid function 
Vi will be referred to as A-function (or V-function) if there exist grid nodes m 
and n such that m < n; v„i > fm-i and Vi > Vi^i (w,„ < Vm-i and w,; < Vi-i 
for the V-function) if i < m; Vi — const if m < i < n; w„ > Vn+i and Vi > u^+i 
(u„ < Vn+i and Vi < Vi+i for the V-function) if i > n. Simply stated, the 
A-function (or V-function) is a scalar grid function Vi that has only one gen- 
eralized local maximum |50) (or generalized local minimum |50| . respectively). 
The set of /^-functions (or ?7-functions) is the union of A-functions (V-functions, 
respectively) and the set of monotone functions. 

Hereinafter, for the sake of convenience, a monotonicity preserving scheme 
will be referred to as G- monotone for short, a scheme monotone in terms of 
Harten et al. [5S] will be referred to as H-monotone, a scheme will be referred 
to as S-monotone if a maximum principle holds for this scheme (see, e.g., [55] . 
[551 p. 183], [TD], [H]), and a scheme being monotone from the standpoint of 
Ostapenko [50] will be referred to as GO-monotone [lO) . 

For studying G-monotonicity of non-linear schemes the notion of total vari- 
ation (TV, see, e.g., [H], [36], [40]) turns out to be an useful tool, since any 
total variation diminishing (TVD) scheme is G-monotone [211 P- 168], [26], [40l 
p. 110]. H-monotone as well as TVD schemes appear to be attractive for at 
least the following reasons: (i) H-monotone schemes are TVD [26] and, hence, 
G-monotone [26]; (ii) any H-monotone scheme, being TVD, converges to the 
physically relevant solution, i.e. solutions of H-monotone schemes do satisfy the 
entropy condition [5S], [5S]; (iii) The notion of a TVD method is sufficient to 
prove convergence [101 p. 148]; (iv) There exist simple sufficient conditions for 
a scalar scheme to be TVD [HI p. 169]. Besides, it is widely believed that TVD 
methods are free from spurious oscillations (e.g., [13], [26], [32], [40], [41], [48]); 
and thus, for the above-stated reasons, TVD schemes are in common practice 
(see, e.g., [21], [36], [37], [40], [44], m, M, [SZ], M, M)- 

Let us note that Harten's theorem relative to G-monotonicity of H-monotone 
schemes (i.e., H- monotonicity =^ TVD G-monotonicity) was proven in [261 p. 
360] for a specific class of schemes approximating a 1-D scalar partial differential 
equation (PDE), and hence it does not always happen that an H-monotone 
scheme is G-monotone. Actually, let us consider the following linear scheme: 

Vi = aiVi^i + P^Vi + ^iVi+i, (4) 

where at = 1 ~ 2e, f3^ ~ J^ = e a.t i ~ 1, 2, 3, . . ., and a, = /3j = e, 7^ = 1 — 2e 
at i = 0, —1, —2, . . ., (0 < e < 0.25). Scheme ^ is H-monotone, since a^, 
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7i > Vi. Considering the monotone increasing function Vi, namely Vi = for 
all i < and Vi — 1 for all j > 0, we obtain that Vi = for all? < 0, = 1 — 2e, 
vi — 2e, and Vi = 1 for all i > 1. We note that the grid function Vi will not be 
monotone. Thus, H-monotonicity is not sufficient as well as not necessary (see, 
e.g., [ini Example 4.2]) condition for a scheme to be free of spurious oscillations. 

Furthermore, as demonstrated in [lOj . a conservative scheme, H- monotone 
(hence, consistent with the entropy condition [25]), TVD (hence, G- monotone 
[55]), and S-monotone, be it non-linear [TUl pp. 1578-1580] or even hnear with 
constant coefficients [TUl p. 1561], can produce spurious oscillations comparable 
with the size of the jump-discontinuity (cf. [IH])- Thus, the notions of TVD, 
S-monotonicity, and G-monotonicity are not sufficient for a numerical scheme to 
be non-oscillatory. As is shown in [10] , the notion of GO-monotonicity is a very 
helpful tool for the construction of non-oscillatory schemes. However, a GO- 
monotone scheme can be not TVD and, hence, not S-monotone, since the 
notion TVD can be viewed as a concept within S-monotonicity [10) . Hence, if a 
numerical scheme will be GO-monotone as well as S-monotone (GOS-monotone 
for short), then this scheme will be stable and, in general, free of spurious 
oscillations [lO] . 

An approach to investigate non-linear difference schemes for S-monotonicity 
in terms of corresponding variational schemes was suggested in [10], [H]. The 
advantage of such an approach is that the variational scheme will always be 
linear and, hence, enables the investigation of the monotonicity for nonlinear 
operators using linear patterns. It is proven for the case of explicit schemes 
that the S-monotonicity of a variational scheme will guarantee that its original 
scheme will also be S-monotone [TOl . Analogous theorem for the case of implicit 
schemes can be found in Section 12.11 Theorem [2] Since a variational scheme 
carries such an important properties of its original scheme as S-monotonicity 
and GO-monotonicity ([10], see also Section [2.11 Theorem [2]), the following 
definition will be useful in the investigation on scheme monotonicity. 

Definition 1 A numerical scheme is termed variationally monotone if its vari- 
ational scheme is monotone. 

The notion of TV turns out to be an effective tool, [40l p. 148], for studying 
the stability of non-linear schemes. Actually, the following property 

||AA(v + ,5v) -A/'(v)ll < (l + aA<)l|^v|| (5) 

is sufficient for stability of a two-step method [40 , however it is, in general, 
difficult to obtain. Here At denotes the time increment, a is a constant inde- 
pendent of At as Ai ^ 0, V and 6v are any two grid functions {Sv will often 
be referred to as the variation of the grid function v), TV denotes the scheme 
operator. At the same time, the stability of linearized version of the non-linear 
scheme is generally not sufficient to prove convergence [17], [ID]- Instead, the 
TV-stability adopted in [27] (see also [JO] s. 8.3.5]) makes it possible to prove 
convergence (to say, TV-convergence) of non-linear scalar schemes with ease. 
However, the TVD property is a purely scalar notion that cannot, in general. 
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be extended for non-linear systems of equations, as the true solution itself is 
usually not TVD [H], gO]. Moreover, one can see in [TOl pp. 1578-1581] that a 
TVD scheme can be non-convergent in, at least, Loo, in spite that the scheme 
is TV-stable (see, e.g., gOl Theorem 12.2]). Such a phenomenon is caused by 
the fact that TV is not a norm, but a semi-norm. 

Nowadays, there exists a few methods for stability analysis of some classes 
of nonlinear difference schemes approximating systems of PDEs (see, e.g., [18], 
PU] . [10], mi, [H], [55] and references therein). It is noted in [ID] that the 
problem of stability analysis is still one of the most burning problems, because 
of the absence of its complete solution. In particular, as noted in [TH] in this 
connection, the vast majority of difference schemes, currently in use, have still 
not been analyzed. LeVeque [40] noted as well that, in general, no numerical 
method for non-linear systems of equations has been proven to be stable. There 
is not even a proof that the first-order Godunov method converges on general 
systems of non-linear conservation laws [40] p. 340]. Thus, a different approach 
to testing scheme stability must be adopted to prove convergence of non-linear 
schemes for systems of PDEs. The notion of scheme in variations (or variational 
scheme [10], [11]) has, in all likelihood, much potential to be an effective tool 
for studying stability of nonlinear schemes. Such an approach goes back to 
the one suggested by Lyapunov (1892), namely, to investigate stability by the 
first approximation. This idea has long been exploited for investigation of the 
stability of motion [T^]) as well as the stability of difference equations pp] . 
We establish the notion that the stability of a scheme in variations implies the 
stability of its original scheme (see Section 12.11 Theorem [3] and Remark [4]) . 

Aiming to demonstrate potentialities such notions as 'variational scheme' 
and 'GOS-monotonicity' in construction numerical schemes, we will develop a 
novel version of the central Lax-Friedrichs (LxF) scheme (e.g., [ST] p. 170]) with 
second-order accuracy in space and time. The stability of this scheme is proven 
in Section [4?2l We restrict our proof mainly to the case when the solutions to 
(|T]) are smooth. Notice, such a property as smoothness is peculiar not only to 
vanishing- viscosity (e.g., [H], [JD]) solutions of a hyperbolic system. Since the 
stiffness parameter, r, in the system ([1]) is like a viscosity jlD], a solution to the 
system of conservation laws (jlj is expected to be smooth provided that input 
data are sufficiently smooth. This specific type of systems is interesting itself, 
as such systems are not uncommon in practice. For a detailed discussion on this 
subject, including the sufficient conditions for the global existence of smooth 
solutions, see [24] . 

An extensive literature is devoted to central schemes, since these schemes 
are attractive for various reasons: no Riemann solvers, characteristic decompo- 
sitions, complicated flux splittings, etc., must be involved in construction of a 
central scheme (see, e.g., [B], [37], [35], [3D], [5T], [S3| and references therein), 
and hence such schemes can be implemented as a black-box solvers for general 
systems of conservation laws |37| . Let us, however, note that the numerical 
domain of dependence [40l p. 69] for a central scheme approximating, e.g., a 
scalar transport equation coincides with the numerical domain of dependence for 
a standard explicit scheme approximating diffusion equations |40' p. 67]. Such 
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a property is inherent to central schemes in contrast to, e.g., the first-order 
upwind schemes (40l p. 73]. Hence, central schemes do not satisfy the long 
known principle (e.g., [31 p. 304]) that derivatives must be correctly treated 
using type-dependent differences, and hence there is a risk for every central 
scheme to exhibit spurious solutions. The results of simulations in [JS] can be 
seen as an illustration of the last assertion. Notice, all versions of the, so called, 
Nessyahu-Tadmor (NT) central scheme, in spite of sufficiently small CFL num- 
ber (Cr — 0.475), exhibit spurious oscillations in contrast to the second-order 
upwind scheme {Cr = 0.95). The first order LxF scheme exhibits the excessive 
numerical viscosity. Thus, the central scheme should be chosen with great care 
to reflect the true solution and to avoid significant but spurious peculiarities in 
numerical solutions. 

Let us note that LxF scheme - the forerunner for central schemes [6], [37] 
- does not produce spurious oscillations. While, from the pioneering works of 
Nessyahu and Tadmor [48] and on, the higher order versions of LxF scheme can 
produce spurious oscillations. The reason has to do with a negative numerical 
viscosity introduced to obtain a higher order accurate scheme. Let us illustrate 
it with the scheme (135) in [9], which is 0{At + (Aa;)^) accurate. We rewrite 
this scheme to read 

i vi^f -viVo.5 f(vivi)-f(vr) _ 

2 0.25Ai Ax 

Ax^ vf - 2vr^o.5 + v-Vi Ax^ dl,, - 
At Aa;2 AAt Ax ' ^ ' 

where df denotes the derivative of the interpolant at x = Xi. Notice, the 
second term in the right-hand side of ([5]) is, in fact, the negative numerical 
viscosity. Without this term. Scheme (O would be of the first order, namely 
LxF scheme. Let us note that there is a possibility to increase the scheme's 
order of accuracy, up to 0{{At)'^ + (Ax)^), by introducing into Scheme ([6|) an 
additional non-negative numerical viscosity. Such an approach is similar to the 
vanishing viscosity method [21], [40], and hence possesses its advantages, yet 
it appears to be free of the disadvantages of this method, since the additional 
viscosity term is not artificial. With this approach, the second order scheme is 
developed in Section 14.21 where sufficient conditions for stability as well as a 
necessary condition for S-monotonicity of the scheme are found. The developed 
scheme is tested on several conservation laws in Section \5\ 

A stable numerical scheme may yield spurious results when applied to a stiff 
hyperbolic system with relaxation (see, e.g., [2], [5], [7], [H], [14], [30], [54] . 
[55]). Specifically, spurious numerical solution phenomena may occur when un- 
derresolved numerical schemes (i.e., insufficient spatial and temporal resolution) 
are used (e.g., [2], [30], [32], [46]). However, during a computation, the stiff- 
ness parameter may be very small, and, hence, to resolve the small stiffness 
parameter, we need a huge number of time and spatial increments, making the 
computation impractical. Hence, we are interested to solve the system, ([1]), with 
underresolved numerical schemes. It is significant that for relaxation systems a 
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numerical scheme must possess a discrete analogy to the continuous asymptotic 
limit, because any scheme violating the correct asymptotic limit leads to spu- 
rious or poor solutions (see, e.g., [H], [SU], [3T], [35], [S]). Most methods for 
solving such systems can be described as operator splitting ones, [Tl], or meth- 
ods of fractional steps, [7]. After operator splitting, one solves the advection 
homogeneous system, and then the ordinary differential equations associated 
with the source terms. As reported in [23], this approach is well suited for the 
stiff systems. Let us however note that the schemes based on the operator- 
splitting techniques are of the first order in time, excluding rare cases such as, 
e.g., the Strang splitting [3D]. Thus, following Pareschi [5TJ p. 1396], we con- 
clude that such schemes are, in general, not robust. Fortunately, as applied 
to System ([!]), the second order schemes can be constructed on the basis of 
operator-splitting techniques with ease. We are mainly concerned with such an 
approach in Section 14.31 



2 Monotonicity and stability of difference schemes 

2.1 Non-linear schemes 

We consider a nonlinear implicit scheme 

H.(yi,...,yAf) =x„ zecj = {l,2,...,M}, (7) 

where y^GL = denotes the sought-for vector- valued function of grid nodes, 
XiGL = denotes the prescribed vector- valued function of grid nodes, Hi= 
{Hi,i, . . . , Hi,N}^ is a vector-valued function with the range belonging to . 
If we introduce the additional notation y={yf , . . . , yj^f}'^ , x={x^, . . . , xjj^'^ , 
H={Hf , . . . , nlj} , then the scheme (O can be represented in the form 

H(y)-x, xeL*^yei^'^. (8) 

Theorem 2 Let a nonlinear implicit scheme ^ be written in the form 
where x. Eilx C i^^, fix denotes a closed and bounded convex set. Let the 
mapping H in (0) have a strong Frechet derivative (strong F-derivative 14 9\ 
item 3.2.9]), H'(y), at every y S int{ily) provided that H (ilj,) = fix, and 
let H' (y) be nonsingular. Then, for any x, x+(5x G ^Ix the scheme will be 
S-monotone if its variational difference scheme is S-monotone. 

Proof. The scheme (jS]) can be seen, [TT] p. 1126], as a two- node implicit 
scheme, and its variational scheme becomes 

<5x = H' (y) . <5y, H'(y)^^^. (9) 

As H' (y) is nonsingular, we may rewrite © in the form 

^y ^ (H')"' (y) ■ <5x. (10) 
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Then, by ^ 



\\Sy\\^ (H')-'(y)-<5x < (H')-'(y) l|<5x|| 



Let ^ be S-monotone, i.e. let 

ll'^yll < ll^xii . 

In view of (HI]) and JTID we obtain [11, p. 1126] that 



(H')"' (y) 



< 1, Vy eQy C L 



M 



In view of the inverse function theorem ■49], item 5.2.1] we obtain from 



-M 



y = H-i(x), X er!, C L*^ y G^iy C L 
By virtue of the mean- value theorem item 3.2.3] we obtain from ([H)) 



(11) 
(12) 

(13) 

that 
(14) 



IMyll 



|H-i (x + dx) -H-i (x)|| < sup (H-i)'(x + idx) 

0<t<l 



Idxil 



< sup 



H 



(H(y)) 



lldxll 



In view of the inverse function theorem [49', item 5.2.1] we can write 

(H-i)'(H(y)) = (H')-'(y). 
By virtue of and we obtain from that 

IMyll < IMxil . 

The inequality (fT7|) manifests the prove of the theorem. ■ 
Theorem 3 Let a non-linear explicit scheme be written in the form 
v = H(v), V e L, V e c L, 



(15) 

(16) 
(17) 



(18) 



where fl denotes a closed and bounded convex set in a linear vector space L. 
Then for any v, v + i5v S $1 the scheme will be stable if its variational scheme 
is stable. 



Proof. The variational scheme corresponding to the scheme reads 

aH(v) 



(5v = H'(v) • 5v, H'(v) 



(19) 



The hnear scheme (HH) wiU be stable gUl p. 145] if ||H'(v)l| < 1 + aAt for all 
V e 57, that is 

sup||H'(v)|| < 1 + aAi. (20) 



v60 



By virtue of the mean-value theorem [151 item 3.2.3] we obtain from 

||H(v + 5v) -H(v)|| < sup ||H'(v + CH|| Pv|| < sup||H'(v)|| ||Jv||. (21) 

o<c<i veo 

In view of (PO)) we conclude from (PT|) that the inequality ([5]) for will be 
fulfilled, and hence the original non-linear scheme will be stable. ■ 

Remark 4 r/ieorem[5| can be reformulated for implicit schemes with ease. The 
proof of this theorem is identical to the proof of Theorem\^ 

2.2 Linear schemes 

We will consider explicit linear schemes on a uniform grid with time step At 
and spatial mesh size Ax. In view of the CFL condition [ID], we assume for 
the explicit schemes, that Ai — 0(Ax). Moreover, we will also assume that 
Ax = O (At), since a central scheme generates a conditional approximation to 
Eq. Il]) (see Section [3l). In such a case, the following inequahties will be valid, 
for sufficiently small At and Ax, 

:^oAt < Ax < /ipAt, i^o,fJ.Q ~ const, < < ^1q. (22) 

Notice, for hyperbolic problems it is often assumed that At and Ax are related 
in a fixed manner (e.g., [101 p- 140], [571 P- 120]), i.e. it is assumed that At 
and Ax fulfill a more strong condition than (P^ . 

Let V = {. . . , vf , . . .}"^ (or A = {Aij}) be a partitioned [43j vector (or ma- 
trix, respectively), then we shall denote by (v) the ordinary vector obtained from 
V (or by (A) the ordinary matrix obtained from A, respectively) by removing 
its partitions. It is easy to see that 

||v|U^max||v,:|U = ||(v)||^. (23) 

To start with, we obtain necessary conditions for some class of linear schemes 
to be GOS-monotone. We consider the following explicit homogeneous scheme 

z» = X! ' ' 2» ' ^ (24) 

j 

where L denotes the linear vector space with the orthonormal basis {b;}^, bi 
= {1,0,..., 0}^, b2 = {0, 1, . . . , 0}^, bM = {0,0,..., if; B„- = {B^j } is a 
square matrix. It is assumed that any constant (i.e., Zi — yj = c = const) is a 
solution to Then, in view of (gH), we find that 

J^B,, =1, yi (25) 

3 

will be the necessary conditions for ((24l) to be G-monotone (cf. [lOl p. 1560]). 
Here and in what follows, I denotes the identity operator. 
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Theorem 5 Let an explicit linear scheme he written in the form \24^ , and let 
any constant he a solution to {24^ . If {24-^ is GOS-monotone, then the diagonal 
elements, B'^^ , of the matrices By = | Bf^ } are non-negative, i.e. 

B^^>0, yi,j,k, (26) 

and 

B^^ is a ^ — function of i for all j and k. (27) 

Proof. We consider ([M]) when y-j — yjhi, I = 1,2, ...,M, where yj is a scalar 

value. Let {z{ ^, z\ ^, .... Zjv/.i}^ be the left-hand side of ((24)) under = yjhi. 
Then we obtain from ((24l) the following system of decoupled scalar equalities 



= k,l^ 1,2,..., M. (28) 



In view of Corollary 3.14 in [TU], the scheme will be S-monotone iff 

^|B^'|<1, yi,k,l. (29) 



In view of (1^51) we have 



and hence 



Y.B'y^^l, Vz,fc, (30) 



Y,\B^^\>^, ^hk. (31) 

j 

By virtue of pTjl and using under k = I we obtain that 

^|B,f;^-| = l, Vz,A:. (32) 

i 

Thus, (|30p and (15^ must be valid simultaneously. It is possible iff all coefficients 
B^j" comply with ((26)) . 



To prove ((27|) we consider (|28|) under k = I. Let m be the scheme matrix 
column number and Smj denote the Kronecker delta. Assuming that the scheme 
will be GO-monotone, it will transform yj = Smj into a /i-function as Smj is 
a A-function of j. Then we obtain from ((^5]) . under k — I, that z^'^ = -Bf^. 
Hence, Bf^ is a ^-function of i, \/k,m. ■ 

Consider the special case of the scheme ((24|) . namely, in ((24)) depends 
on a square matrix 

By = (Ai) , Vi, j, (33) 

where A,; is similar [43', p. 119] to a diagonal matrix A^, i.e. there exists a 
non-singular matrix Si such that 

• A, . S, = A,=diag[xlxl . . . , Af } . (34) 
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It is assumed that B^j = if j ^ Ji = {j : i — kL<j<i + kj^}, where fc^,, fc^ 
= const > 0. Notice, it is not assumed here that any constant (i.e., Zi = = 
c = const) is bound to be a solution to ([M)) . The foUowing notation is used: 

y = {. . . ,yj, . . .}^ , yj = SJ^ ■ yj, y ={..., yj, .. .}^ , y^^^ = S^^ . y^-, 

~ _ / — T — T —T —T "1^ 

\ • . • ■ y.i^i — fe^ _i J yi i — k^ I ■ . • J Yi^i+kn 1 yi,i+fcn + l ' ' ' ' / ' 

B={Bij}, By = Sj "'^ • By • Si, Bi= |. . . , Bij_i, By , By+i, . . .} . (35) 
The stabihty for ([M)) provided can be addressed by the foUowing. 

Lemma 6 Consider an explicit scheme that can be written in the form \24^ un- 
der hSi^). \S4-^ . Let Si = s (Aj) be the spectrum of A^, and ip^j (A) be represented 
by an absolutely convergent power series at each point A G s^. Let By = in 
^24^ if j ^ Ji = {j ■ i — kL<j<i + kn}. Then the scheme will be stable if 

max^ kij (^)| < 1> Vi, (36) 

j 

||(Sri - StI) . Sj^ < e = const, Vi, Vj G Ji. (37) 
Proof. It is easy to see that 

Sri^{l+(S-i-S-i).S,}.S-i. (38) 
Then, in view of (j37|) . we find Vi, Vj G Ji 

||Sr^.y,L<||l+(Sri-S7i).S,||^||STi.y,.||^< 

{1 + II (Sri - STi) . S, 11^} ||STi . y,.||^ < (1 + e) IIS7I . y, 11^ . (39) 
By virtue of ([55)) . we rewrite ([M]) to read 

z, = Sri . z, - ^ By • (Sri • y^) = B, • y, = (B,;) • (y,) , Vz, (40) 
3 

where (Bi) and (y^) denote the ordinary matrix and vector obtained from Bj 
and yi, respectively, by removing the partitions. Notice, By = if j ^ Ji, since 
By = for these j. Thus, it can be written that By — S,^^ • By • if j ^ Ji, 
and hence y j ^ — yj j ~ SJ^ ■ yj = yj (Vj ^ Ji). In view of ((40|l we obtain that 

||z,|U^||S,ri.z.||^<||(B.)||^||(y.)IU, Vz, (41) 

The norm IKyi)!!^^ in (PT|) can be estimated by virtue of and 

II (y.) IL = llydL < (l + ©) max \\Sj' ■ y, 11^ = (1 + e) llylL , V*. (42) 
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Let us estimate ||(Bi)||^ in ([iTjl . In view of = S ■ ^ • • S^. It can be 

verified, by induction with respect to n, that (A^)" = S^^ • (A.^)" • S^. Then, in 
view of Theorem 11.2.2 and Theorem 11.2.4 in [13], we find 

B,, ^ Sri . B,, . S. = (p,, (S-i • A, • S,) = (A,) . (43) 
Thus, Bij can be written in the form 

B,,^dtag{Al^,Al,...,Afl}, A^- = (a)') , fc = l,2,...,M. (44) 
In view of (jH]), we find that 

By virtue of dH]), gS]), and ^ we obtain from (gT]) that 

||Sri.z,||^<(H-e)||y||^, Vi. (46) 

Since 

l|z|loo = mfx||Sr^ -z^ll^, (47) 
we obtain, in view of (|46|) . ([47|. that 

l|zlL^Noo<(i + 0)llylloo = (i + 0)llyL- (48) 

The last inequahty estabhshes Lemma HI ■ 

We consider the following explicit linear scheme 

vr+i^^Br^-vJ, n>0, (49) 
j 

where 

B^-j"-!;"^^' , V.,,n, (50) 

-^i = : i ~ kL < j < i + kn} , kL.kji^ const, Vi, n, (51) 

fci, /ci^, denote the non-negative integer constants being independent of t, a;. 
Ax, and At. It is assumed that the matrix-valued function A — A{x,t) is 
Lipschitz-continuous and A is diagonizable, i.e. for A"= A (xi, tn) there exists 
a non-singular matrix S" ~ S (xi,t„) such that 

(Sr)-' • Ar • Sr = Ar^d*ag {Ar\ A^^ . . . , Ar^^"} , Vz, n. (52) 

Let us note that even if B^- = iy9^ (A") in ([50)1 and Lemma [5] be valid for 
the linear scheme (|1^ with = (Ai) at every time step, the scheme (|1^ will 
be "locally stable" only, i.e. any growth in error is, at most, order O (At) in 
one time step. However, we cannot, in general, show on the basis of that 

||v^^||^^ <Ct||v°||^, Ct = const, (53) 
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where ||-||^^ and ||-||^ denote some norms, v" = |. . . , (v")"^ , ■ • - j , Nt denotes 
the time level corresponding to time T — NxAt over which we wish to compute. 
The reason is that the vector norm in (|48|) depends on the time level i„, and 
hence we maynot apply (|48|) recursively to obtain ((53|) . The stability of the 
system can be addressed by the following. 

Theorem 7 Consider an explicit scheme that can be written in the form 
under (5d^-(5^, where the functions iffj (A) and A{x,t) are both Lipschitz- 
continuous. Let there exist Axq > such that the function (A) in h5U\) 
can be represented by an absolutely convergent power series at each point of 
the spectrum s" = s (A") \li,n, Vj G Ji, VAx < Axq, and let the matrix-valued 
functions S {x, t), and S^^ {x, t) in I5i?|) can be taken such that the matrix-valued 

■ S" will be 



functions 



(Sf)-' - (S")"' {x) ■ S" (x) and (S,)-' (t) - (Sf) 



Lipschitz-continuous in space and, respectively, time Vi, n. Let 



< f3_i = const, 



|S"| 



< Pq = const, \/j,n. 



Then the scheme 114 9\ j will be stable, i.e. \53\) will be valid, if 
max^ (A)| < 1, yi,n. 



Proof. Let B" — ipf, (A"), and let us rewrite (I49p to read 



Vi, n, 



where 



j 

j 

First, let us estimate the norm, hn (•), of v" : 

hn (V") 



I v'^ 1 1 = max 

I M oo 



Since \{Sfy^ - (S")"^ -S" {x) and [(S^)"^ {t) - (S^')"^ 
continuous in space and time, respectively, we may write 



(54) 
(55) 

(56) 
(57) 
(58) 

(59) 

■ S" are Lipschitz- 



■ s 



i+l 

■s 



< PiAx, (3i=const, \/i,n, (60) 
< (32At, (32— const, \/i,n. (61) 



By virtue of ([60|) . we find 



on 



< P^Ax, ^3 = const, yi,n, Vj £ J,;, (62) 
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where — /J^ max(fci, fc^). We assume for the expUcit scheme (UHl), that 
Aa; = 0{At) (i.e. 3 Aio > 0, 3 ao > such that Aa; < aoAt VAt < Ato). 
Then we find by virtue of ([5^ that 



(63) 



The inequahty (j63p coincides with the assumption p7p in Lemma [B] under O 
/34At. Then, in view of LemmajHl we obtain for the scheme ([ST]) , that 



(SD" 



< 



[l + /34Ai]max (S^"' • = [1 + /34At] /i„ (v") . 



(64) 



Let us now estimate the norm /i„(v"). Since (/?^- (A), A(x,t) are both 
Lipschitz continuous, we may write 

ll^r, (A,") - (ADL ^ «i II - A,"||^ , Vi,n, V^- e J., (65) 

||A^ - A^ll^ < aa -Xil < asAx, Vi,n, VjgJ^, (66) 

where as = aa max(fci, fc^), ai, aa = const. By virtue of ([65)) . ([66]) . and 
assuming that Ax = O (At) , we obtain 

||B^^-B^^||^^||^^^(A7)-^^^(Ar)||^<a4At, V*,n, VjeJ., (67) 

where 0:4 — aoaias — const. We obtain from (|58p that 

Whence, by virtue of (p7)) and we obtain 



(68) 



< asAimax 



(S;') '-v^ = a5At/i„ (v") , Vz,n, (69) 



where as = /?_4/?oa4 max(fcL, fc^j) — const. By virtue of ([5S)) . ([M]) . and 
we obtain 

hn (v"+i) < [1 + /3Ai] /i„ (v") , /? = /34 + as = consi, Vn. (70) 
It is easy to see that 

{s-+y ^ {i + {{s^y - is-y') ■ sr} • (sr)"' , (71) 

whence, by virtue of ([6T|) . we find 



/i„+i (v"+i) = max (S'/+i) '-v; 



< 
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max 



\ / oo 



< 



(1 + At) max 



In view of ^ and ((TlD, we find 

K+i (v"+i) < [1 + lAtf K (v") , 7 = max (a, P^) , Vn. 
Applying ((73|) recursively gives 

hNr (v^^-) < (1 + jAtf''^ ho (v°) < Ct/^o (v°) , Cr - cxp {2^T) . 



(72) 
(73) 
(74) 



The inequalities in (l74|) prove the theorem. ■ 

Let us consider the case when the operator B,^- in (I49p depends on a matrix 
A" belonging to a set of pairwise commutative diagonizable matrices: 



B" - (A") , A" . = A^ • A", Vz, j, n, k, m. 



(75) 



In such a case, the S-monotonicity of the system (|49|. can be addressed by the 
following. 

Theorem 8 Consider an explicit scheme that can be written in the form 
provided ^75^ . Let ip^j (A) in (TSp cott, 6e represented by an absolutely convergent 
power series at each point of the spectrum = s (-^j ) Ji '^^ T'^en the scheme 
|^ff[ ) wi/Z 6e S-monotone iff 



AGs" ' ' 



(76) 



Proof. As A" belongs to the set of pair-wise permutable diagonizable ma- 
trices, the matrices of the set are simultaneously similar to diagonal matrices 
|431 p. 318], i.e., there exists a non-singular matrix S such that 



S-i-A«.S = A, 



'J=diag {aJ'\ A"'^, . . . , A"'^^| , Vj, i 



(77) 



where A"'™ denotes the m-th eigenvalue of A". The following notation is used: 



v" = S"^ • v" B 



■ Br 



S, B" = {b:;}, B"^{B^^}. (78) 
By virtue of ([75)) . we rewrite (|^ to read 

-„+i^S-i.vri=^B:;.v^. (79) 



Using v" = |. . . , (v^)'^ , ■ • - j , we rewrite ([79| to read 



v"+i = B"-v", 



v"+l^{...,(v^l)^...}' 



(80) 
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In view of ([50]) we obtain that 

/i(v"+i) = ||(v"+i)||^ < 



B 



(81) 



where j and (v") denote the ordinary matrix and vector obtained from 

B" and v", respectively, by removing the partitions. Let us estimate the norm 

of ^B"^ in HHID. Since (A) can be represented by an absolutely convergent 

power series at each point A G s" — s (A") , we find, in view of Theorem 11.2.2 
and Theorem 11.2.4 in 43J, that 



B 



B^ . S = (S- 



" (A") 



Thus, Bi,- can be written in the form 



(82) 



B,,=rf^ag{A:;/^A:^^...,A5^}, a:^'^ = ^» (a;^'^) , fc = l,2,...,M. (83) 
In view of (1551). we obtain that 



B' 



max I max A"/ 

7 \ fe=l,...,J\/ ^ I 



Whence, in view of (I76p . we find 



B' 



< 1, Vn. 



(84) 



(85) 



The vector norm h (•) in (I81|l docs not depend on time level. Then, in view of 
Proposition 3.2 in llj, the inequality ([85]) proves Theorem [8] ■ 



Proposition 9 // {24-^ is a variational scheme, then k25\) is, in general, not 
valid. Notice, Lemma and Theorem are proven without assumption i25\) . 
However, in addition to the Lipschitz- continuity ofA(x,t) (see 133\) and i50\) ). 
it is assumed in Lemma\^ and Theorem^ that some functions of S (x,t) (see 
{^P; (E3i ) oltr also Lipschitz-continuous. Let us note that the stability of a 
linear scheme can often be proven without assumption I125\) as well as without 
additional assumptions on the continuity. To demonstrate it, let us generalize 
the theorem of Fnednchs (1954) (see, e.g., p. 120], fBO, p. 374]) to be 
applicable to variational schemes. We consider the following difference scheme 



y"+i (x) = J2 Bfe (^) • y" ix + kAx) , y" G B^t G 



cMxM 



(86) 



k——m 



where x G (— oo, oo), B" is a symmetric and non-negative matrix. It is assumed 
that y" (x) and BJJ [x] are periodic (with the period equal to 1) functions of x. 
Let 



= y" (x + J Ax) , B^,^. =B]^{x + jAx) , 



(87) 
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and let 



(u, v) = y [u^ (x) -v (x)] dx, l|u|| = ^/M. 



If there exist ci, C2 = const such that 



C2 



2m + 1 



Az, 



(88) 



(89) 



then the scheme is stable. Notice, it is not assumed that (x) = I. 

The proof is very little different from the proof when B^ (x) — I. Actu- 
ally, in view of I186\) and the first inequality in I189\) . we obtain 



,n+l I 



(B^ . yl?, y"+i) < 0.5 ^ [(B^ ■ y^, y^) + {B^ ■ y"+\ y"+i)] 
fc fc 

< 0.5 J2 (Bfc • Yfc , Yfc ) + 0.5 (1 + ci Ax) ||y"+i f . (90) 



Since y" (x) and BJ! (x) are periodic functions, we obtain, by virtue of (2^ 
fgP|) . and that 



(1 - „,c,At) ||y"+i f < J2 (Bfc,fc • yfc , yfc) + E ((B^ - ^^I'c) • yfc , Xfc) < 



^ (B^ . y", y") + ^-£^ Ax ^ (y^ y,") = (1 + Mo (ci + C2) At) ||y"f . (91) 



follows from !i91\} that 



„+l||2 1 + Mo (ci + C2) At 



1 - /ioCiAt 



(92) 



Let Ato = const such that l—fiQCiAto > 0. In particular, let Ato = 0.5/ (/ipCi). 
Then, for all At < Ato the following inequality will be valid 



,n+l I 



< (1 + C3 A<) ||y" f , C3 = 2^0 (2ci + C2) = const. (93) 



The inequality in i93\) proves Proposition\^ 

Notice, in practice, we often deal with systems for which the matrices B^ in 
Scheme are not symmetric. Let us generahze the theorem of Friedrichs, [S71 
p. 120], [6OI p. 374], to be apphcable to variational schemes with non-symmetric 
matrices B^. 

It is assumed that there exist matrix- valued functions W {x, t) and (x, t) 
such that W = W^, W • = (W • A^)^ , 

Ao(u,u) < (W-u,u) < Ao(u,u), Ao,Ao = const, < Ao < Ao, (94) 
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||W"+' - W"|| < aoAi, ao^ const, W"=W(a:,i„), (95) 
||A^-B^||w„ <^-^Ax, /3o = const, A^-Afe(x,t„), (96) 

II - A^^^- 1|^„ < 6oAx, 60 = const, A'^^ = A^ (x + jAa;, i„) , (97) 

where ||A|| is the norm of an operator A on a Hilbert space equipped by the 
scalar product ([55)1 . while ||A||^ is the norm of an operator A on a Hilbert 
space equipped by the foUowing scalar product 



(u,v)w = (W-u,v), ||u||w = ^(u,u)w. (98) 

The generalization can be addressed by the following theorem, where the 
notation jHT]), dHH]), and jMl) wiU be used. 



Theorem 10 Consider an explicit scheme that can he written in the form ISf^) . 
where y" G R^'^, BJ! G R^^^"*^^ are periodic (with the period equal to 1) functions 
of X. Let W(a;,i) he symmetric, i.e. W = W"^, and positive (W > in 
terms of ^94^ ) matrix-valued function, and let i96]) . ^97\ l he valid for the case of 
symmetrizable, i.e. W"-A^ = (W"-A^')'^, and periodic matrix-valued function 
A^(^)- //EfcAfc(x,t) - I anrf 

Ai (u, u) < ((W • Afe) -u, u) < Ai (u, u) , Ai, Ai = const, < Ai < Ai, (99) 

then Scheme i86\) will he stahle 



Proof. Multiplying (|86p by W", and by y"'"''"'^, we obtain, after integrating 
both sides, that 



^ llw 



< J2 (AI! •y^,y"+i)^„ ((BI! - A^) .y^,y"+i)^„ . (100) 



Since W" • A^' is symmetric and, in view of positive, and since J2k AJJ — I, 
we write: 

^ (A^ • y^, y"+i)^„ < 0.5 ^ (A^ ■ y^, y^)^„ + 0.5 ||y"+i ||^„ . (101) 
fc fc 

By virtue of ([W)) . periodicity of and y^?, and since J^k Afe = I, we obtain: 

E (Afc ■ y^yl?)w" = E i^h ■ y^y^)w"+E ((^1^ - ^h) ■ ylNy^Ow- ^ 



E (A^ • y", y")w. + boAx E (yl?, y^)w.. < (1 + biAx) ||y"|l^,. , (102) 

k k 



where bi — (2m + 1) 6o = const. In view of ()96|) and the assumption of period- 
icity, we obtain for the last term in the right-hand side of (jlOOp that 

E ((B^ - K) ■ yl, y"+i)^„ < E IKBfe - A^) ■ ylU. ||y"+i ||^,. < 
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[3,^x lly" ||^„ ||y"+i ||^„ < 0.5/3oA:r |ly" ||^,. + 0.5/?oAx ||y"+i . (103) 
After elementary transformations we find from (|100p - (|103p that 



l|y""^lll.< '^^^"^y^" lly"ll^^- (104) 



Let Axq = const such that 1 — /JqAxq > 0. In particular, let Axq = 0.5/ /3q. 
Then, in view of (|22p . for all Aa; < Axq the following inequality will be valid 

||y"+iw" < (l + 62At)||y"||^„, 62 = (2/3o + fei) - const. (105) 
By virtue of and ([Ml), we find that 

||y"+iw"+i = (W"+i •y"+\y"+i) - (W".y"+i,y«+i) + 
((W"+i-W") -y^+Sy^+i) < ||y"+i||^„+«oAt(y"+\y"+i) < 

(l + 63At)||y"+i||5^,., 63 = ^ = const. (106) 
Then, in view of (|105p and (|106p . we write 

||y"+lw"+i ^ (l + c^*)l|y"llw- max (62, 63) = const (107) 

The inequality in (I107P proves Theorem [101 ■ 

3 Monotone piecewise cubics in construction 
of central schemes 

In this section wo consider some theoretical aspects for high-order interpolation 
and employment of monotone piecewise cubics (e.g., |16| . [35]) in construc- 
tion of monotone central schemes. By virtue of the operator-splitting idea (see 
also LOS in [SHj), the following chain of equations corresponds to the problem 

Q 

2^ + ^f(u) = 0, t„<t<t„+o.5, u (2;, <„) = u" (x) , (108) 

1 dw. 1_ 

--^ = -q.{n), t„+o.5 < i < t„+i, u(x,t„+o.5) = u"+°-5 (x) . (109) 

2 at T 

Using the central differencing, we write 



du 

9t 



n-l-0.25 n 



f-f x-r 0.25At 

1 — 171 + 0.125, X — Xi-\-0,5 



dx 



O ({Mf) , (110) 



fn-l-0.125 _ fn-l-0.125 

--^ h 
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By virtue of (|110p - (|llip we approximate (|108p on the cell [x^, x^+i] x i„+o.25] 
by the following difference equation 



'gr+7-"'-gr"-"'), (112) 



n+0.25 _ /'„n+0.125 „Ti+0.125 



where v"^'^, g"^'' are the grid functions. As usually, the mathematical treat- 
ment for the second step (i.e., on the cell [xi_o.5 7 2:^+0.5] x [tn+o.25, ^n+o.s]) of a 
staggered scheme will, in general, not be included in the text, because it is quite 
similar to the one for the first step. 

Considering that (|112p approximate pOSp with the accuracy 0((Ax)^ + 
(Ai)^), the next problem is to approximate vf_^Q 5 and g"^"'^^^ in such a way 
as to retain the accuracy of the approximation. For instance, the following 
approximations 



'i+O.S 



0.5(vr+vr+i)+o((Aa;)'), gr+°-^'' - f (v^ + O (At) , (113) 



leads to the staggered form of the famed LxF scheme that is of the first-order 
approximation (see, e.g., [HI p. 170]). One way to obtain a higher-order scheme 
is to use a higher order interpolation. At the same time it is required of the 
interpolant to be monotonicity preserving. Notice, the classic cubic spline does 
not possess such a property (see Figure [1^). Let us consider the problem of 
high-order interpolation of v"_^q g in (|112p with closer inspection 

Let p = p (a;) = {p^ (x) , . . . (x) , . . . ,p'^ {x)}^ be a component-wise 
monotone piecewise cubic interpolant (e.g., (TB], [SS]), and let 



p., = p (xj) , p- = p' (xi) , Ap,, = Pj+1 - Pj, 

p,=A..— , p,^,=B..— , (114) 

where p^ denotes the derivative of the interpolant at a; = Xi. The diagonal 
matrices and in (|114p are defined as follows 

A,^dtag{alal...,aT}, = dtag {p], PI . . . , PT} ■ (US) 

The cubic interpolant, p — p{x), is component-wise monotone on [xi,Xi+i] iff 
one of the following conditions (e.g., [16], [35]) is satisfied: 



a? - 1)^ -f' (a? - 1 



) (pt - 1) + (p'; -if --i («■■ + /3- - 2) < 0, (116) 

a? + /3- < 3, a^, > 0, /3f > 0, Vi, k. (117) 

As reported in |35| . the necessary and sufficient conditions for monotonicity of 
a piecewise cubic interpolant originally given by Ferguson and Miller (1969), 
and independently, by Fritsch and Carlson [16^. The region of monotonicity 
is shown in Figure [T]d. The results of implementing a monotone piecewise 
cubic interpolation when compared with the classic cubic spline interpolation. 
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Figure 1: Monotone piecewise cubic interpolation, (a) Interpolation of a 1-D 
tabulated function. Circles: prescribed tabulated values; Dashed line: classic 
cubic spline; Solid line: monotone piecewise cubic, (b) Necessary and suffi- 
cient conditions for monotonicity. Horizontal hatching: region of monotonicity; 
Unshaded: cubic is non-monotone. 



are depicted in Figure [T^. We note (Figure [T^) that the constructed function 
produces monotone interpolation and this function coincides with the classic 
cubic spline at some sections where the classic cubic spline is monotone. 

Using the cubic segment of the piecewise cubic interpolant, p = p(a;), 
(see, e.g., |16j . |35j ) for x G [x^jXi+i], we obtain the following interpolation 
formula 

p,+o.5 = 0.5 (p, + p,+i) - ^ (p^+i - pO + O {{^xY) . (118) 

If p [x) has a continuous fourth derivative, then r = 4 in (|118p . see e.g. [331 P- 
111]. However, the exact value of p^ in (jllSp is, in general, unknown, and hence 
to construct numerical schemes, employing formulae similar to (jllSp . the value 
of derivatives p^ must be estimated. 

Using ()118|) and the second formula in (I113P we obtain from (I112p the fol- 
lowing scheme 

<t^f = 0.5 (vr + - ^ (dr_,, - ar) - ^ '^"H^'^"'^ (119) 

where d" denotes the derivative of the interpolant at x — Xi. In view of (jllSp 
and the second formula in (|113p . the local truncation error [JDl p. 142], ip, on a 
sufficiently smooth solution u{x,t) to (|108p is found to be 

= O (At) + O (^^) + O ((At)2 -I- (Axf) . (120) 

In view of (|120p we conclude that the scheme (|119p generates a conditional 
approximation, because it approximates (|108p only if {AxY /At — > as Ax 
and At ^ 0. Let d" be approximated with the accuracy O {{AxY), then 
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the value of r in (|120p can be calculated (see Section [HI Proposition [TT|) by the 
following formula 

r = min(4, s+1). (121) 

Interestingly, since (|119p provides the conditional approximation, the order of 
accuracy depends on the pathway taken by Aa; and At as Ax — > and Ai 0. 
Actually, there exists a pathway such that At is proportional to (Ax)'^ and 
the CFL condition is fulfilled provided /i > 1 and Aa: < Axq, where Axq is a 
positive value. If we take n = 1 and s > 1, then we obtain from (|120p that 
the scheme (|119p is of the first-order. If = 2 and s > 3, then (|119p is of the 
second-order. However, if /i = 2 and s = 2, then, in view of (|120p and (|121[) . the 
scheme ()119p is of the first-order. Moreover, under fi = 2 and s — 2, the scheme 
will be of the first-order even if g"'^°'^^^ in (I112|) will be approximated with the 
accuracy 0{{Aty). It seems likely that Example 6 in [37] can be seen as an 
illustration of the last assertion. The Nessyahu-Tadmor (NT) scheme with the 
second-order approximation of d" is used [37] to solve a Burgers-type equation. 
Since At = 0((Aa;) ) [37], the NT scheme is of the first-order, and hence it 
can be the main reason for the scheme to exhibit the smeared discontinuity 
computed in [371 Fig. 6.22]. 

The approximation of derivatives can be done by the following three 
steps [12]: (i) an initialization of the derivatives pj; (ii) the choice of subregion 
of monotonicity; (iii) modification of the initialized derivatives p[ to produce a 
monotone interpolant. 

The matter of initialization of the derivatives is the most subtle issue of this 
algorithm. Actually, the approximation of pj must, in general, be done with 
accuracy 0((Ax)^) to obtain the second-order scheme when At is proportional 
to (Ax) , inasmuch as central schemes generate a conditional approximation. 
Thus, using the two-point or the three-point (centered) difference formula (e.g. 
[35] . [51]) we obtain, in general, the first-order scheme. The so called limiter 
functions [35j lead, in general, to a low-order scheme as these limiters are often 
0(Aa:;) or 0((Ax)^) accurate. Performing the initialization of the derivatives pj 
in the interpolation formula (|118p by the classic cubic spline interpolation [56j . 
we obtain the approximation, which is 0((Aa::)'^) accurate (e.g., [3^, [35]), and 
hence, in general, the second-order scheme. The same accuracy, 0{{Ax)^), can 
be achieved by using the four-point approximation |35j . However, the efficiency 
of the algorithm based on the classic cubic spline interpolation is comparable 
with the one based on the four-point approximation, as the number of multi- 
plications and divisions (as well as additions and subtractions) per one node is 
approximately the same for both algorithms. We will use the classic cubic spline 
interpolation for the initialization of the derivatives P'^ in the interpolation for- 
mula (jllSp . as it is based on the tridiagonal algorithm, which is 'the rare case 
of an algorithm that, in practice, is more robust than theory says it should be' 

m- 

Obviously, for each interval [a:;.i,Xi+i] in which the initialized derivatives p^, 
Pi^i such that at least one point (a*^, /?*^) does not belong to the region of 
monotonicity (|116p - (|117p . the derivatives p,^, p'^^i must be modified to p'^, p^^i 
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such that the point {a^ , f3^ ) will be in the region of monotonicity. The modifi- 
cation of the initialized derivatives, would be much simplified if we take a square 
as a subregion of monotonicity. In connection with this, we will make use the 
subregions of monotonicity represented in the following form: 

< a,^ < 4K, < /3f < 4K, k, (122) 

where H is a monotonicity parameter. Obviously, the condition (|122l) is sufficient 
for the monotonicity (see Figure [Jd) provided that < H < 0.75. 

Let us now find necessary and sufficient conditions for (|118p to be G-monotone. 
By virtue of (|114p . the interpolation formula (|118|) can be rewritten to read 

Pi+0.5 = (0.51+ I ■p^ + (0.51 I -p.+i. (123) 

The coefficients of (|123p will be non-negative ijf |/3j — ai\ < 4. Hence (|118p 
will be G-monotone iff (|122p will be valid provided < K < 1. Notice, there 
is no any contradiction between the sufficient conditions, (|122l) provided < 
H < 0.75, for the interpolant, p = p(a;), to be monotone through the interval 
[ijjXi+i], and the necessary and sufficient conditions, (|122p provided < H 
< 1, for the scheme (I123|) to be G-monotone. In the latter case the interpolant, 
p = p (x), may, in general, be non-monotone, however at the point i + 0.5 the 
value of an arbitrary component of Pi+0.5 will be between the corresponding 
components of Pi and p^+i. 

To fulfill the conditions of monotonicity (|122p , the modification of derivatives 
p- = {pf ,pf , . . . can be done by the following algorithm suggested, in 

fact, by Fritsch and Carlson [T^] (see also |35j): 

~ 4Kminmod(Ati, A*^), minmod(pf , 5*^), H = const, (124) 

where A^ = {Pi^i — p^) XAx, the function minmod(a;, y) is defined (e.g., [55] , 
[37], [15], [5T], [S]) as follows 

minmod(2:, y) = — [sgn{x) + sgn{y)] min (|a::| , |y|) . (125) 

4 Construction of first- and second-order cen- 
tral schemes 

Central difference schemes with first- and second-order accuracy are introduced 
in this section. The construction of the central schemes is based on: (i) Vari- 
ational GOS-monotonicity notion, (ii) Monotone piecewise cubic interpolation 
(e.g., [H], [35]), (in) Operator-splitting techniques (see also LOS in [59]). 

4.1 First-order central schemes 

Let us note that instead of point values, v"_)_q 5, employed in the construction of 
the scheme (jll2p . it can be used the cell averages (e.g., [5], [37], [ID]) calculated 
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on the basis of the monotone piecewise cubics. In such a case we obtain, 
instead of (|118p . the following interpolation formula 

Ax 

Pj+0.5 = 0.5 (p, + p,;+i) - K— (p'^+i - , (126) 
where >c — 2/'3. The region of monotonicity in this case will also be 

< A, < 4HI, < < 4KI, < H < 1, Vi. (127) 

Notice, the interpolation formula (|126p coincides with pisp under x — 1. Thus, 
in view of the interpolation formula (|126p . the staggered scheme (I112p is written 
to read 



V 



n+0.25 



,«,s -0Mv?«+v:.)-.^(d»„-df)-^I<^4-IWl, (128) 

where d" denotes the derivative of the interpolant at x — Xi, the range of values 
for the parameter is the segment < ><• < 1. As usually, the mathematical 
treatments for the second step of the staggered scheme will, in general, not be 
included in the text, because the second step is quite similar to (|128p . li k =1 
(or >f = 0), then Scheme (|128p coincides with the scheme (|119p (or with the LxF 
scheme, respectively). As it was shown above, the scheme (|128p is of the first 
order provided At = O (Ax). In such a case, since the source terms can be, in 
general, stiff (i.e., r ^ 1), it is natural to use the following first-order implicit 
scheme for ()109p . 

The first order central scheme (|128p - (ll29p based on operator-splitting techniques 
will be abbreviated to as COSl. 

Let us investigate the stability of Scheme (I128p . It is assumed that the 
vector-valued function u = u (x, t) is Gateaux- (or Frechet-) differentiable on 
the convex set ilxt C R x [0, +oo), and its derivative is bounded on Q.xt- It 
is also assumed that the operator A (= (9f (u) /du) is Frechet-differentiable 
on the convex set fiu C M}'^ , and its derivative is bounded on fiu- Hence A 
= A (x, t) will be Gateaux- (or, respectively, Frechet-) differentiable [49l p. 62] 
and its derivative will be bounded on ^Ixt^ and hence A (x, t) will be Lipschitz- 
continuous on [49l p. 70]. Since the Jacobian matrix A = A{x,t) in ([2]) 
possesses M linearly independent eigenvectors. A" (= A [xi, t„)) is similar to a 
diagonal matrix [43j . i.e. there exists a non-singular matrix S" = S {xi, t„) such 
that 

(S^r' ■A^-S^^diag{K'\K'\...X'^'}, Vz,n. (130) 

The right and left eigenvectors of Af~ A (xi, i„) can be defined in such a way 
that Sf = S{xi,tn) will be the matrix having the right eigenvectors as its 
columns, and the rows of (S")^^ = {xi.tn) will be the left eigenvectors 
[39l p. 62]. For Theorem [7] to be used, it must be proven that (S")~^ — 
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(S")"^ (a;)] • S" (x) and [(S,)"^ (t) - (Sf • will be Lipschitz-continuous 

in space and, respectively, time Vi, n. Since each of the above functions depends 
on one parameter (x or t) only, it can be done with ease for strictly hyperbolic 
systems, i.e. when the eigenvalues of the operator A (= di (u) /du) in ([2]) 
will be all distinct [401 P- 2]. In this case the proof follows from the long- 
known results of perturbation theory for simple eigenvalues [62^ p. 67] (see also 
Theorem 2.1 in 4 ). For a detailed discussion on the derivatives (sensitivities) of 
eigenvectors of matrix-valued functions depending on several parameters when 
the eigenvalues are multiple see [4], |63j . 

By virtue of (|114p , the second term in right-hand side of (|128p can be written 
in the form 



(divi - dr) = f (Br 



Then, the variational scheme corresponding to ()128p is the following 



(131) 



r n+0.25 



o.5(<5vr+5vr+o 

At 



r • {5^ - ^vJVO + (Ar • <5vr - Ar,i • 5vIVi 



(132) 



2Ax 

where = diag {-D"i, -D"2) • ■ ■ > m] = B" - By virtue of (fTTT)) . we find 
that -4KI < < 4HI, and hence -SHI < JUf < 8KI. Thus, we may write 
that 



IIJBrils <8H. 



(133) 



Considering that v" in |l3 



is Lipschitz-continuous on ^xt, we write 

Ci, = consi. (134) 



It is assumed, see (|22l) . that there exists = const such that Ax < aoAt for 
a sufficiently small At. Then, by virtue of (|133p and (|134p . and since < >ir, H 
< 1, we find the following estimation for the second term in right-hand side of 



- 8 II » 



'i+l|l2 



||<5I 



< aoCyAt. (135) 



In view of (|135p . the scheme (|132p will be stable if the following scheme will be 
stable (e.g., [Si pp. 390-392]) 



'i+0.5 



0.5 ( I+^l 



Ax 



0.5 (1 + EH • <5v; 
• (5v" + 0.5 



0.5 (I - E^Vi) • '^v^Vi 



4 ' Ax 



'i+i- 



We write, in view of Q and (fT27)) . that 

At 



s(E«)c 



Aa 



At 



Vi, n. 



(136) 



(137) 
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Hence, by virtue of Theorem [7] we find that the scheme (|136p will be stable if 
max 0.5(|1 + A| + |1-A|) ^ 1, \/i,n. (138) 

Aes(E7) 

We obtain from p38p the following condition for the stability of the variational 
scheme (|136p : 

>^K + a^i, a = ^^^, (139) 

where Cr denotes the Courant number. Thus, in view of Theorem[3]the scheme 
will be stable if will be valid. 

Notice, (fT39)) was obtained on the basis that Ef = f + ^jA^ is diagoniz- 
able. Such an assumption should be verified for a concrete problem. It might be 
well to point out that this assumption has a rigorous basis if A" is diagonizable 
and is a scalar matrix, i.e. = 46*"!, 6*" = const (-H < 6*" < H). In such 
a case E" will be diagonizable. Thus, for instance, LxF scheme will be stable 
if Cr ^ 1. In the case, when D" is not a scalar matrix, whereas the Jacobian 
matrix A in (l2|) is symmetric, the condition (|139p for stability of (|128l) can be 
found with ease by virtue of Proposition [9l In a more general case, when E" 
is not diagonizable as well as A is not symmetric, the stability of (|128|) can be 
investigated by virtue of Theorem 1101 

Let us find a necessary condition for the variational S-monotonicity (see 
Definition[l]) of the COSl scheme, p28)) - (fT29)) . Considering (fT32)) on a constant 
(v" = C = const, yi,n), we obtain 

^vr++o°5'' - 0.5 (<5vr + 5v,IVi) + |b • {S^r- S^^^,) + 

Here, in (|140p . Ac is a diagonizable matrix such that its spectrum s (Ac) C 
[— Amax, Amax], scc ([2]). In vicw of (|13ip . D in (|140p can be taken at will, as 
v" = C = const, vf_(_j — vf = and d" = d"_^;^ = 0. Aiming to obtain the 
necessary condition for (|132p to be S-monotone, it is assumed that D = is a 
scalar matrix with C G [— 4K,4K]. In view of Theorem [51 the scheme (jl40l) will 
be S-monotone ijf 



max 

|C|<4K, AGs(Ae 



1 kC At ^ 

- + — H —A 

2 8 2Ax 



1 kC At ^ 

2 ^ ~8~ ^ 2A^ 



^ 1. (141) 



By virtue of (|14ip . we find that (|139p is the necessary condition for the varia- 
tional S-monotonicity of (|128p . 

The variational scheme corresponding to (|129p reads 

5vf+i = 6^^+°-' + (vf+i) • 5vf+i, (142) 

where G (V) = dq (V) /dV. Let us rewrite (IT42|) to read 

<5vr+i = (l-— G(vr+i)l '.Jv^+o-^ (143) 
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Let G will be a normal matrix. In such a case, in view of fTU^, Theorem 3.3] the 
variational scheme (|143p will be S-monotone if 



max^ 

k 



1 - ^4 (vr+^) 



< 1, Vi,n, (144) 



where denotes the k — th eigenvalue of the matrix G. In view of ([3]), the 
inequality (|144p is valid, and hence the variational scheme (I142p will be uncon- 
ditionally S-monotone. 

Thus, the COSl scheme, (|128p - (|129p . will be variationally S-monotone only 
if pM]) will be valid. 

4.2 Second-order central scheme 

In this section, the second-order scheme for (|108p is developed by approximating 
v^+o 5 and g"+°-^^^ in pT^ with the accuracy 0{{Axf + (Atf). The sufficient 
conditions for stability as well as the necessary condition for S-monotonicity of 
this scheme are found. 

Using Taylor series expansion, we write 



At 



dt 

By virtue of the PDE system, pOS]) . we find 

di di du „ di di dn „ / 9f \ ^ dn 



O (At^) . (145) 



dt du dt du du dx \ du ) dx ^ 

Using the interpolation formula (|126p and the formulae (|145p - (|146p . we obtain 
from (|112p the following second order central scheme 

vr/o°f = 0.5 (vr+1 + vD - (dr+i - dr) + 



(At) 



2 



8Ax 

Since d" is the derivative of the interpolant a,t x = Xi, the third term in the 
right-hand side of (|147p can be seen as the non-negative numerical viscosity 
introduced into the first order scheme (|128p . Owing to this term, the scheme 
(|147p is 0{{Ax) + (At) ) accurate. Since the source terms in Q can be, in 
general, stiff (i.e., r ^ 1), it is natural to use the following second-order implicit 
Runge-Kutta scheme for (|109p . since this scheme possesses a discrete analogy 
to the continuous asymptotic limit, 

vr+i=vr+°-^+7qK+°-^5), 7^^- (148) 
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Scheme (|147p - (|148p . based on operator-splitting techniques, wih be abbreviated 
to as C0S2. 

Let us find the sufficient conditions for stability of Scheme (|147p . It is 
assumed that the following inequalities are valid in a suitable norm 



9vr 



< Pa Pv" II , ttA, /3a = const. 



(149) 



In view of l|114p , the variational scheme corresponding to (|147l) is the following 



r n+0.25 



= o.5(<5vr+5vr+o+? (vr-vr+o'^-^ror +>r-('^vr-<5vr+o 



(Ml 



{K(A,IVi)' 



8Ax2 



(5v^Vi-'5vr) + 



2Ax 



(Ar • - A^+i • 5vr+o 



(150) 



In view of Q, p33)) . p34l) . and pl9)l . Scheme ^50]) wiU be stable if the following 
scheme will be stable. 



(Ml 

8Aa;2 



(A-Vi) 

At 
2Aa; 

We rewrite (|15ip to read 



- (AD' 



(Ar • (5vr - A^, . 5v^0 



^vr++o°5'' - 0.5 (I + ED • + 0.5 (I - E^+O . <5vr+i, 



where 



4 



E 



4 (Ax) 
(At)' 



(ADi)' 



(AD' -A, 



At 

Aa; 



A", 



(ADi)'-B,-(AD' 



An 
J A 

4"' 4(Ax)' L^"-'^^^ -'\^ ^+1' 

We write, in view of © and (fT77)) . that s (ED C [-A^, A^], where 



Ab = >fH — 



/AtA„ 



At 

H + -r— Ainax, Vz, U. 

Ax 



(151) 
(152) 

(153) 
(154) 

(155) 



V Ax 

Hence, by virtue of Theorem [7] we find that the scheme (|152p will be stable if 

AtA 



(>r-c2)H + a < 1, a = 



Ax 



(156) 
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Thus, in view of Theorem [31 the scheme p47p will be stable if (|156p will be 
valid. 

By analogy with the COSl scheme, we find the necessary conditions for the 
variational S-monotonicity (see Definition [T|) of the C0S2 scheme, (I147p - (|148p . 
assuming that v" = C = const. The variational scheme corresponding to p47p 
under pTi)l . is the following 



e «+0.25 



= 0.5 + 5vr+i) + > • (5vr - <5v,'Vi) 



8 

At 
2Ax 



(157) 



8 {AxY 

where = A(C), A (u) = df (u) /du, D = (I- In view of Theorem El the 
scheme (|157p will be S-monotone Zj(f 



max 

|C|<4N, Aes(Ae) 



8 



3 (Ax)' 



At 
2 Ax' 



8 



8 (Ax)' 



At . 
2 Ax' 



^ 1. 



(158) 



By virtue of (|158[) we find that (|156p will be the necessary condition for the 
variational S-monotonicity (see Definition [ij of the non-linear scheme (|147p . 

The variational scheme corresponding to (|148l) (under v" = C = const) 
reads 



Ge - G (C) , 



(159) 



where G (u) — dq (u) /9u. Let G will be a normal matrix. In such a case, 
since all eigenvalues of Gc have non-positive real parts, see ([3]), the first step 
of the scheme (|148p will be S-monotone (see [lOl Theorem 3.3], Theorem [2]). It 
remains to prove that the scheme (|159p . taken as a whole, will be S-monotone 
under the same condition. Eliminating ^v""*""'^^ we obtain that 



I-?G. 



I 7Gc 

In view of [TUl Theorem 3.3], the scheme (|160p will be S-monotone if 

1 



max 

k 



l-74(C)(l-i^fc(C)) 



< 1, VC en^, 



(160) 



(161) 



where (C) denotes the fc-th eigenvalue of the matrix Gc. Since all eigen- 
values of Gc have non-positive real parts, the variational scheme (|160p will be 
unconditionally S-monotone. 

Thus, the C0S2 scheme, (|147p - (|148p . will be variationally S-monotone only 
if pro]) will be valid. 



29 



4.3 Second order schemes based on operator splitting tech- 
niques 

Using the first order schemes (|128l) and (I129|) , it can be constructed (see Section 
[6l Proposition [12]) a scheme approximating ^ with the accuracy 0{{Ax) + 
{At) ). ActuaUy, smce 



/A-i-\2 / q2 \ n+0.25 / a J-^2 / o2 \ "+0-25 

(At) ( d u\ [At) (O U 



32 V 9t^ J t+Q.5 32 V dt^ 
the scheme will be as follows 



0{Ax{At)^), (162) 



vr°-^^=vr + ^q(vr+"-^^), (163) 



(d^+°-25 - dr+0-25) _ f^iAIl±L_J_Ll!:^ L_ (164) 



Hence, the stability as well as monotonicity behavior of (|164p and (|163p can 
be interesting itself. As demonstrated in Section [5l the behavior of Scheme 
(|164|) is similar to the one of the NT scheme, i.e. it can produce spurious 
oscillations if CFL number is not sufficiently small. Let us develop another 
scheme approximating H]) with the accuracy 0{{Ax)'^ + (At)'^) and such that its 
components (after operator splitting) will be of the second order. It can be done 
on the basis of the second order scheme (|163l) - (|164p with ease. Actually, adding 
to and subtracting from Equation (|229p (see SectionlHl PropositionfT^ the same 
quantity l|162p . we obtain (after operator splitting) the following scheme, instead 
of (fTCD-(fTMD. 

"'^'+2r'^' 32 U^V. ' ^ ' 

{Atf fd^Y^"'' A^ fK+r^)-f(vr°"^) 

32 U^V.+o.s 2 Ax ■ ^ ) 

Thus, Scheme ()165p as well as Scheme (|166p are of the second order, and Scheme 
(|165p - (|166p . taken as a whole, is of the second order as well. 

Using Taylor series expansion, and central differencing, we find 

\ n+0.125 

-7^'-'' - -7 + f [^) +0{{At)'). (168) 
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Considering the equation in (jlOSp at i„+o.25 < t < tn+0.5, we obtain, in view of 
(fT46l) . that 



d^u _ d fd{\ _ d fdf 



(169) 



where A —di /du. Then 



dx \ dx 



1 

Ax 



-| n+0.25 
i+0.5 



fA"+o-25)'.dr 



o 



By virtue of pUS)) . P^ - p7D)) . we rewrite Scheme p^ - pM]) to read 



(170) 



^n+0.125 _ ^ri+0.25 



— ('n"+°-12'^ + n+0.25\ 



,n+0.25 
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n+0.125 



= 0.5 (v: 



'i+0.5 



n+0.25 
i+l 



.^n+0.25\ 



(171) 



8Ax 



f An+0.25\^ jn+0.25 



(a: 



ri+0.25\ 



iri+0.25 



Atf(vl!+^ 



T!+0.25\ 



n+0.25\ 



"aT^^ (^^2) 

Notice, Scheme (|172l) coincides, in fact, with Scheme (|147l) . however there is 
a difference between (|17ip and (|148p . In spite of the fact that both of these 
imphcit Runge-Kutta schemes, (|17ip and (|148[) . are of the second order, only 
the above combination, i.e. Scheme (|17ip - (|172p . is of the second order as a 
whole. Thus, even if a higher order implicit Runge-Kutta scheme will be used 
to approximate (I109p . nevertheless, there is a risk that Equation ^ will be 
approximated with the first order in time. 

It can be proven, by analogy with the proof in Section 14.21 that p56p is 
the sufficient condition for stability as well as the necessary condition for S- 
monotonicity of Scheme (|17ip - (|172p . 

There exists one more problem associated with the stiffness, i.e. r ^ 1, 
of (|109p . To demonstrate it, let us consider the following system of ordinary 
differential equations (ODEs) with a relaxation operator in the right-hand side: 



dX (t) _ ip {X, Y) 

dt ^ T 

dYjt) ^ ip{X,Y) 

dt T 



{X-Y), X(0)-Xo, 



{Y-x), y(o) = ro. 



(173) 
(174) 
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The system will be considered at the interval < i < At. The following notation 
will be used: = X (I'At). If (p{X,Y) = 1, then the approximation, Xi, to 
Xi, in the case when the first order scheme (|129p is used, will be 

X*,=Zo + ^f^, (175) 

where Zq = 0.5{Xo + Yq), A = 2At/T. In the case of the second order scheme, 
(|148[) . the approximation to Xi will be 

xr^Z,+ (176) 
^ 1 + A + 0.5A^ 

The analytic solution, Xi, will be the following 

X, = Z, + (177) 
exp (A) 

Notice, in the case of underresolved numerical scheme we obtain that A ^ 1. 
A comparison between (I177P and (|176p shows that there is a need to develop an 
implicit Runge-Kutta scheme of higher order accuracy, and such that the scheme 
(based on the operator splitting techniques), taken as a whole, will be, at least, 
of the second order. However, the higher order will be the implicit Runge-Kutta 
scheme, the more time-consuming procedure will be obtained. To tide over this 
problem, a semi-analytic approach could be useful. Let us demonstrate one, as 
applied to System (fT73)) - (fT74)) . Let 0^ = ip{X^,Y^). Integrating (fT73)) -p74 l) . 
and using the midpoint rule, we obtain 



Xi — Xq — — - 

T 



Yi-Yo = - — 

T 



At 

0.5 



J {X - Y) dt + O (^{Atf'^ , (178) 





At 



J {Y -X)dt + {{Atf^ . (179) 







To resolve (|178p - (ll79p . a hnearized version of System (|173p - (|174p can be used, 
namely, instead of (p{X,Y) in (|173p - (|174p . it can be taken 60.5 — const. The 
analytic solution will be the following 

— Zf) Yc] — Zn 

^1 = ^0 + (on A.y ^ ' ^1=^0 + ° ° . . ■ 180 

exp (2&o,5At/T) exp {2&o,5At/T) 

It remains to estimate the value of 6'o.5. Since 60.5 = 0.5{tp {Xq, Yo) + ip {Xi,Yi)) 
+ 0((Ai)^), we obtain the following, in general, non-linear equation in ^0.5 

00.5 = 0.5((p (Xo, Yo) + if (Xi, n)), (181) 

where Xi and Yi are defined via (|180p . 
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5 Exemplification and discussion 



In this section we are mainly concerned with verification of the second order cen- 
tral scheme, (|147|) . To demonstrate the superiority of this scheme over Scheme 
(|128|) being 0{{Ax) + At) accurate, we will use the inviscid Burgers equa- 
tion. The C0S2 scheme, (|147l) - (|148p , is verified using Pember's rarefaction test 
problem [55j . Euler equations of gas dynamics, namely Sod's as well as Lax 
problems, are also used for the verification of the second order central scheme, 

urn- 

5.1 Scalar non-linear equation 

As the first stage in the verification, we will focus on the following scalar version 
of the problem ([T]): 



du d 
dt dx 



f{u) = Q, xGR, 0<t<r,,ax; u{x,t)\^^^ = u^{x). (182) 



To test the first order scheme, (|128p . we will solve the inviscid Burgers equation 
(i.e. / {u) = u'^/2) with the following initial condition 

M(a;,0) = |^°' ^x^^ih^j^^) ' > Hl, ^ const ^ Q. (183) 
The exact solution to ()182|) . ()183|) is given by 

where T — 2S / uq, S = Hr — Hl, 



''^-mq, Hl < X < b, b — UQt + Hl 



b-IiL 

ui{x,t) = <^ uq, b < X < O.duQt + Hr , (185) 

0, X < Hl or X > O.buot + Hr 

( 2S{x-hL) u ^ ^ T 

u,(xA = \ 7^^""' ^^<^^^ , (186) 
0, X < or X > L 



L = 2^S^ + O.buoS {t-T) + Hl- (187) 

The numerical solutions were computed on a uniform grid with spatial in- 
crements of Ax = 0.01, the velocity uq = 1 in p83p . = 0.2, Hr — 1, the 
monotonicity parameter K = 0.5, the Courant number Cr = uoAty Ax = 0.5, 
and the parameter x — 1 in (|128p . The results of simulation are depicted with 
the exact solution in Figure [21 We note (Figure [2|) that the first order scheme, 
(|128p . exhibits a typical second-order nature, however spurious solutions are 
produced by the scheme. To obtain necessary conditions for the variational 
GOS-monotonicity (see Definition [T]) of the COSl scheme, (|128p . we consider 



33 





t=0 / 




t=l 


1=2 






0.5 





Figure 2: Inviscid Burgers equation. The scheme COSl {>c = 1) versus the 
analytical solution. Crosses: numerical solution; Solid line: analytical solution 
and initial data. C,. = H = 0.5, Aa; = 0.01. 

the scheme in variations (|132[) . corresponding to (I128p . on a constant (w" = C 
= const, Vi,n), i.e. (|140p . By virtue of Theorem [8] and Theorem [5l we find 
that, in the case of the Burgers equation, Scheme (|140|) will be GOS-monotone 
only if (|139p will be fulfilled. Since all of the coefficients of (|140p will be non- 
negative (the scalar version is considered) under (|139p . the scheme in question, 
(|128p . will be H- monotone (and hence TVD and G-monotone) only if (|139p will 
be valid. Notice, the numerical simulations were performed with such values of 
the parameter >c, Courant number, Cr, and monotonicity parameter, H, that 
(|139p was not violated. As it can be seen in Figure [H the boundary maximum 
principle is not violated by the scheme, i.e., the maximum positive values of the 
dependent variable, v, occur at the boundary t = 0. It is interesting that the 
spurious solution (see Figure [2]) produced by the scheme COSl has the mono- 
tonicity property |26j , since no new local extrema in x are created as well as the 
value of a local minimum is non-decreasing and the value of a local maximum 
is non-increasing. 

Let us note that the problem of building free-of-spurious-oscillations schemes 
is, in general, unsettled up to the present. Even the best modern high-resolution 
schemes can produce spurious oscillations, and these oscillations are often of 
ENO type (see, e.g., [53] and references therein). We found that the oscillations 
produced by the COSl scheme, p28p . are of ENO type, namely their amplitude 
decreases rapidly with decreasing the time-increment At, and the oscillations 
virtually disappear under a relatively low Courant number, Cr < 0.15. However, 
the reduction of the Courant number causes some smearing of the solution. The 
spurious oscillations (see Figure [2]) can be eradicated without reduction Courant 
number, but decreasing the parameter x. Particularly, the spurious oscillations 
disappear ii >c — 2/3, Cr = 0.5, however, this introduces more numerical 
smearing than in the case of the Courant number reduction. Satisfactory results 
are obtained under x — 0.82 {Cr = 0.5). The results of simulations are not 
depicted here. 

To gain insight to why the scheme COSl, ()128p . can exhibit spurious so- 
lutions, let us consider the, so called, first differential approximation of this 



34 



scheme ([T71 p. 45], [BHl P- 376]; see also 'modified equations' in |T7I p. 45], 
[40] . [44]). As reported in [17], [60], this heuristic method was originaUy pre- 
sented by Hirt (1968) (see [TT] p. 45]) as well as by Shokin and Yanenko (1968) 
(see [501 P- 376]), and has since been widely employed in the development of 
stable difference schemes for PDEs. 

We found that the local truncation error, -0, for the scheme COSl can be 
written in the following form 

{l~^)iAxf d^u{x,t) At d'fju) 
^ 4At dx^ 4 dtdx 

O + (Atf + (Axf'^ . (188) 

By virtue of (|188p . we find the first differential approximation of the scheme 
COSl 

du df{u) At d f du{x,ty 



dt dx 4 9a; V dx 

n2 



, (189) 



where B ^ {1 - k) {Ax/ At) - . The term in right-hand side of iflSQ]) wiU 
be dissipative if 



Ax^ ^ 



(1 - >c) j - > 0, =^ Ct <\-^. (190) 

Thus, the scheme COSl, (|128p . is non-dissipative under x = 1, and hence can 
produce spurious oscillations. Notice, ii x — 0.82, then we obtain from ()190p 
that Cr < 0.42. Nevertheless, as it is reported above, satisfactory results can 
be obtained under Cr — 0.5 as well. 

So then, the notion of first differential approximation has enabled us to 
understand that the spurious solutions exhibited by the scheme COSl, (|128p . 
are mainly associated with the negative numerical viscosity introduced to obtain 
the scheme of the second order in space, i.e. 0((Ax)^ + Ai). Let us consider the 
scheme C0S2, (|147|) . approximating (|182p with the accuracy OiiAxf + iAtf). 
Notice, the second order scheme C0S2, ()147p . is nothing more than the scheme 
COSl, (fT28| . with the additional non-negative numerical viscosity. To test the 
scheme C0S2, (|147p . the inviscid Burgers equation was solved under the initial 
condition (|183p . The numerical solutions were computed under the same values 
of parameters as in the case of the scheme COSl, but Cr = 1. The results of 
simulation are depicted with the exact solution in Figured 

We note (Figure [3]) that the scheme C0S2, (I147|) . exhibits a typical second- 
order nature without any spurious oscillations. Increasing the value of H (up 
to H = 1) leads to a minor improvement of the numerical solutions, whereas 
decreasing the value of Cr leads to a mild smearing of the solutions. The 
results of simulations are not depicted here. 
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Figure 3: Inviscid Burgers equation. The scheme C0S2 {k = 1) versus the 
analytical solution. Crosses: numerical solution; Solid line: analytical solution 
and initial data. = 1, H = 0.5, Aa; — 0.01. 



5.2 Hyperbolic conservation laws with relaxation 

Let us consider the model system of hyperbolic conservation laws with relaxation 
developed in ^55^ : 

-^ + -g^az^-Q{w,z), (192) 

where 

Qiw, z) = z — m{u ~ Wo), u ~ w — qaz, (193) 

T denotes the relaxation time of the system, qq, m, a, and uq are constants. The 
Jacobian, A, can be written in the form 

w~qoZ + a -qo{w-qoz) | ^^^^^ 

The system (|19ip - (ll92p has the following frozen 55] characteristic speeds Ai = 
a, X2 = u + a. The equilibrium equation for p9ip - (|192p is 

dw d f 1 
dt dx 

where 



-ui+aw\=Q, (195) 



771 

u*=ii; — 502;*, z^. = [w — ua) . (196) 

1 + mqQ 

The equilibrium characteristic speed A* can be written in the form 

A, {w) - ^tL^ + a. (197) 
1 + mgo 

Pember's rarefaction test problem is to find the solution {w, z} to (|19ip - 
(|192p . and hence the function u — u{x,t), under r ^ 0, and where 

{w,z}^l W,z,(-l)}, x<xo (^gg) 
1 {wr,z^ [wr.)\ , X> Xo 
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< UL = wl - go^* (wl) < Ur = Wr - q^z^ (wr) . (199) 

The analytical solution of this problem can be found in |55| . The parameters of 
the model system are assumed as follows: qo — —1, m = —1, mq = 3, a = ±1, 
T = 10~^. The initial conditions of the rarefaction problem are defined by 

UL = 2, zl ^ m {ul ~ uq) = 1, Wl =ul+ qoZL = 1, (200) 

Ur = 3, ZR = m {ur - Uq) = 0, wr = ur + qoZR = 3. (201) 

The position of the initial discontinuity, xq, is set according to the value of a 
so that the solutions of all the rarefaction problems are identical [53]. Let a 
position, z^, of leading edge or a position, x^, of trailing edge of the rarefaction 
be known (e.g., = 0.85, a;^ = 0.7 in [SS]), then 

xo = - f + a]t^xi~ ( ^ + a) t. (202) 

V 1 + mqa J \l + rnqo ) 

At t = 0.3, under pIDl) - pJT|) we have [S5] 

2, x<0.7 

/8~!o% , 0.7 <2;< 0.85 . (203) 

3, ' a; > 0.85 

The results of simulations, based upon the scheme C0S2, (|147p - (|148l) . under 
different values of the parameter a(a = l, a = — 1) and different values of a 
grid spacing. Ax, are depicted in Figured] 

One can clearly see (Figure 2]) that the scheme C0S2 is free from spurious 
oscillations. Let us also note that the results generated by the scheme C0S2 are 
less accurate in the case of negative value of a than those in the case of positive 
value of a. Specifically, in the numerical solutions produced under a = — 1, 
the representations of the trailing and leading edges of the rarefaction are more 
smeared than those in the solutions produced under a ~ \. Notice, under some 
negative value of a, the frozen and the equilibrium characteristic speeds do not 
all have the same sign. 

5.3 Euler equations of gas dynamics 

In this subsection we apply the second order scheme C0S2, (|147p . to the Euler 
equations of gamma-law gas: 

^"^r^^ +^F(u) = 0, xeR, i>0; u (x, 0) = u" (x) , (204) 

u ={mi,U2,U3}^ = {/9,pt),e}^ , F(u) = +p, (e + p)w}^, (205) 

+ -pt'^, 7 ~ const, (206) 



7-1 2 
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Figure 4: Pember's rarefaction test problem. The second-order scheme C0S2 
{k = 1) versus the analytical solution for u. Dashed line: numerical solution; 
Solid line: analytical solution. Time t = 0.3, Courant number Cr — 1, mono- 
tonicity parameter N = 1. (al): Aa; = 10"^, a = 1; (a2): Ax = 2.5 x 10""', 
a = 1; (bl): Aa; 10"^ a = -1; (b2): Aa; = 2.5 x 10""*, a = -1. 



where p, v, p, e denote the density, velocity, pressure, and total energy respec- 
tively. We consider the Riemann problem subject to Riemann initial data 

u°(a;)^("^ a;<a;o ul,ur = const. (207) 

The analytic solution to the Riemann problem can be found in [40l Sec. 14]. 

Aiming to suppress possible spurious oscillations, the scheme C0S2, (|147|) . 
is modified in such a way as to prevent a violation of the necessary conditions 
(see Theorem [S]) for variational monotonicity of the scheme C0S2, (|147p . The 
modification is to use H ~ H(a;, t) instead of H = const. Leaving out the terms 
equal to zero under v" = C = const in ()147p . we write the reduced variational 
scheme corresponding to the scheme C0S2 in the following form 

Sv^^^C^-Sv^ + 'D2-Sv-^„ (208) 



S^r' = Cl^ol ■ S^T^I + ^7^o:! ■ ^<:o.i, (209) 



where 



cr = ^ (i + I^Ar + bA , Dr = - I^A^+i - Br ) , (210) 



(Ar+,) -Br-iAD -Ar , (211) 
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/-m+0.5 ^ ^ fi , \ n+0.5 , Tan+0.5 \ /r,-, r)\ 



^7^:i = 2 (I - A^^'^+o.i - ■ (213) 
^ [{^to°i)" ■ - {^^:if ■ ■ (214) 

Notice, if the diagonal elements of the matrices C" and D" are non-negative, i.e. 
^n,kk ^ Q jjn,kk ^ then any matrix column of the scheme (|208p consists 
of zeros, but two entries. In such a case this column is a /i-function. Hence, 
at the first stage of modification, it is sufficient to verify Cf'^^ and D,"■'=^ If 
Cf''^'^ < or < for, at least, one value of k, then the value of H at the 

node i (and, if it is necessary, at the neighbor nodes) is reduced up to Kmin > 0. 
Then, the modified values of H, i.e. , are used in the scheme C0S2, (|147p . 
After that we turn to the second stage of modifications. Eliminating Sy'^^^'^ 
from (|208p - (|209p we convert the variational scheme based on staggered spatial 
grid into the following non-staggered form 

^^r^:+' = • + G7 ' + • Sv?^,, (215) 

where 

= c^^c^Li, = ai+S:i-'DU+Tyi^:!-c-, = ixi^i-B-. (216) 

In view of Theorem [Sj we obtain that the foUowing inequalities must be vahd 

i' , , H^- >0, >mm|Fi_^i ,Hil^ j Vi, Vk. (217) 

where p^''^'^ ^ Qn,kk ^ j_j-n,kk jjgjjQ^g ^j^g diagonal elements of the matrices F", 
Gf , and H", respectively. By analogy with the first stage, the value of H" is 
modified at every node, where at least one inequality in (|217p is violated. 

First we solve the shock tube problem (see, e.g., [6], [40], [41]) with Sod's 
initial data: 

f 1 1 f 0.125 1 

(218) 

Following Balaguer and Conde 6J as well as Liu and Tadmor [41] we assume 
that the computational domain is < a; < 1; the point xq is located at the 
middle of the interval [0, 1], i.e. = 0.5; the equations (|204p are integrated up 
to i = 0.16 on a spatial grid with 200 nodes as in [6] and in [41]. The Courant 
number is taken to be Cr = 1 (or Cr = 0.9) in contrast to [B] and [¥r, where the 
simulations were done under At = O.lAx (i.e. 0.13 < Cr < 0.22). The results 
of simulations with the two-stage modification of the monotonicity parameter 
H are depicted in Figure [5l 
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Figure 5: Sod's problem. The scheme C0S2 with two-stage modification under 
Cr = 0.9, H e [0.3, 1] (left column) and = 1, H e [0, 1] (right column) versus 
the analytical solution. Time t = 0.16, spatial increment Aa; = 0.005. 
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The results depicted in Figure [5] are not worse in comparison to the cor- 
responding third-order central results of [411 P- 418] as weU as to the results 
obtained by the fourth-order non-oscillatory scheme in [SI p. 472]. Moreover, 
the scheme C0S2, (|147p . does not produce any spurious oscillations. We con- 
tinue the testing of the modified scheme C0S2, (|147p . using the shock tube 
problem with Lax initial data 
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The results of simulations with the two-stage modification under Cr = 0.9, 
^min = 0.3 (left column) and Hmin = 0.001 (right column) are depicted in 
Figure [H 

One can readily see (Figure O that the scheme C0S2 gives a not worse 
resolution than that obtained by the schemes studied in |4Ij and [6i p. 472] 
without, in fact, spurious oscillations. 

Thus, the algorithm based on Theorem[5]shows a possibility to avoid spurious 
numerical oscillations in a computed solution totally, and hence the second order 
scheme, C0S2, is found to be accurate and robust. However this algorithm can 
lead to a smearing effect that partly reduces accuracy. There are several reasons 
for this side effect of the algorithm. Particularly, the algorithm is not optimal, 
i.e. it gives no more than rough approximations from below to the values of K. 
Moreover, in this algorithm, the value of K at a grid node i is considered as a 
scalar, i.e. it is assumed that the monotonicity parameter, H, is common to all 
coordinates of the vector v^. Thus, the algorithm lacks flexibility, and hence the 
accuracy of the scheme C0S2, (|147p . can be reduced. Notice, the conditions 
of Theorem [S] are not necessary and sufficient conditions for monotonicity (i.e., 
absence of spurious oscillations) of the difference scheme C0S2, and hence the 
total elimination of spurious oscillations on the basis of Theorem [5] can lead to 
an undesirable smearing effect. All this must be given proper weight in designing 
of the difference schemes. 

6 Appendix 

Proposition 11 Let us find the order of accuracy, r, in hll^) if di will be 
approximated by di with the order of accuracy s, i.e. let 

d, = d, + O {{AxY) . (220) 

Let U (x) be sufficiently smooth, then we can write 

U.+i = + C/;+o5^ + lu^'+o, + O ((Ax)') , (221) 

U.. = f/.+05 - C/;+o5^ + ^C/;Vo5 (^)' + O {{Axf) . (222) 
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Figure 6: Lax problem. The scheme C0S2 with two-stage modification under 
Cr = 0.9, K e [0.3, 1] (left column) and K G [0.001, 1] (right column) versus the 
analytical solution. Time t = 0.16, spatial increment Ax = 0.005, Aleph: the 
distribution of K after the second stage of modifications 
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Combining the equalities \221\) and \22Si we obt 



am 



In a similar manner we write: 



4+05 



2 



O ( (Ax)^ ) . (223) 



Aa; 1 ,^„, / Ax^ ^ 



d:+i = UUo5 + UUo5— + ^^^05 [— ) + O ((Aa;)^j , (224) 

d. = - U'U,,^ + ^^1o, + O ((^^)') ■ (225) 

Subtracting the equations \224^ and h225]) . we obtain 



di+i — di 



4+05 



O 



((Ax)^) . (226) 



In view of i226]) and I1220\) we obtain from i223\) the following interpolation 
formula 

U^+o5 = \ (C/4+1 + (di+1 -d)+0 ((Ax)^ + (Aa;)'+') . (227) 

In view of \22'T^ we obtain that r = min (4, s + 1) . 

Proposition 12 As applied to central schemes (see Section \4-.l\ l, we will con- 
struct a second order scheme based on operator- splitting techniques. We will, in 
fact, use the summarized (summed) approximation method f59. Section 9.3] to 
estimate order of approximation. Consider the following equation 

Vu = Viu + T'su = ^ - Lu = 0, VkU=-^- LfcU, fc = 1, 2, (228) 
ot 2 ot 

where Lk is an operator, e.g. a differential operator, a real analytic function, 
etc', acting on u{x,t). We approximate i228\) on the cell [xi,Xi+i] x [tn,in+o.5] 
by the following difference equation with the accuracy 0((Ax) + {At) ) 

jj^ ^ ~ <+0.5 _ ^ n+0.25 _ A,v"+°-25 = 0, (229) 

where it is assumed that the operator L^u is approximated by the operator A^u 
with the accuracy 0{{Ax)'^), i.e. 



ri+0.25 



{L,n)':^°f + o[{Axf). (230) 



In view of the operator splitting idea, to the problem !i229\) there corresponds the 
following chain of difference schemes 

-, ^"+0.25 _ „ 

2 0.25Af ^ ' ^ ' 
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2 0.25At 2 V ; 

One can see from the above that the operator 'PfcU is approximated by IlfcU with 
the accuracy 0{At + (Aa;)^) 

A , / o9 \ n+0.25 

n,u-"f = (7',u)-°f-^(|^)^^^^ +0{iMf + iAxf), (233) 
n2ur_/o"f - (^2u)r+"f + ^ ( . +0 ((At)^ + (Axf) . (234) 



i+0.5 

In view of i233\} - (234^ , the local truncation error \40[ p. 14^], ipj on a suffi- 
ciently smooth solution u{x, t) to i228\) is found to be 

V' = nu = Hiu+Hzu = (235) 

iViu+V2u)'^+lf + O (^{Atf + (Ax)') = O (^{Atf + (Ax)') . (236) 

Thus, implicit Scheme 11231]) together with explicit Scheme i232\) approximate 
i228\} with the second order. 



References 

[1] R. Abgrall and P. L. Roe, High Order Fluctuation Schemes on Triangular 
Meshes, Journal of Scientific Computing, Vol. 19, Nos. 1-3 (2003), pp. 3-36. 

[2] I. Ahmad., M. Berzins, MOL solvers for hyperbolic PDEs with source 
terms, Mathematics and Computers in Simulation 56 (2001) 115-125 

[3] D. A. Anderson, J. C. Tannehill and R. H. Fletcher, Computational Fluid 
Mechanics and Heat Transfer, Hemisphere Fublishing Corporation, New 
York, 1984. 

[4] Alan L. Andrew, K.-W. Eric Chu, Peter Lancaster, Derivatives of eigen- 
values and eigenvectors of matrix functions, SIAM J. Matrix Anal. AppL, 
Vol. 14, No. 4 (1993), pp. 903-926. 

[5] Mark A. Aves, David F. Griffiths, and Desmond J. Higham, Runge-Kutta 
solutions of a hyperbolic conservation law with source term, SIAM J. Sci. 
Comput., Vol. 22, No. 1, pp. 20-38 (2000) 

[6] Angel Balaguer and Carlos Conde, Fourth-order non-oscillatory upwind 
and central schemes for hyperbolic conservation laws, SIAM J. Numer. 
Anal. , Vol. 43, No. 2, pp. 455-473, (2005) 

[7] Frangois Bereux, Lionel Sainsaulieu, A roe-type Riemann solver for hyper- 
bolic systems with relaxation based on time-dependent wave decomposi- 
tion, Numer. Math. 77: 143-185 (1997). 



44 



I. Boglaev, Monotone Iterates for Solving Nonlinear Monotone Difference 
Schemes, Computing 78 (2006), 17-30. 

V. S. Borisov and M Mond, On monotonicity, stability, and construction of 
central schemes for hyperbolic conservation laws with source terms, (2007), 
larXiv:0705. 1109^ ^3 [physics. comp-ph]. 

V. S. Borisov and S. Sorek, - On monotonicity of difference schemes for 
computational physics, SIAM J. Sci. Coniput., Vol. 25, No. 5 (2004), pp. 
1557-1584. 

V. S. Borisov, On discrete maximum principles for linear equation systems 
and monotonicity of difference schemes, SIAM J. Matrix Anal. AppL, Vol. 
24, No. 4 (2003), pp. 1110-1135. 

Russel E. Caflisch, Shi Jin, and Giovanni Russo, Uniformly accurate 
schemes for hyperbolic systems with relaxation, SIAM J. Numer. Anal., 
Vol. 34, No. 1 (1997), pp. 246-281. 

S. R. Chacravarthy, S. Osher, High resolution applications of the Osher 
upwind scheme for the Euler equations, AIAA paper #83-1943, Danver, 
Mass. (19843), pp. 363-372. 

Tao Du, Jing Shi, Zi-Niu Wu, Mixed analytical/numerical method for flow 
equations with a source term. Computers & Fluids 32 (2003) 659-690. 

Duboshin G. N., Theoretical foundations of stability of motion, Moscow 
University, Moscow, 1952 (in Russian). 

F.N. Fritsch and R.E. Carlson, Monotone piecewise cubic interpolation, 
SIAM J. Numer. Anal. 17, No. 2, 238-246 (1980). 

V. G. Ganzha and E. V. Vorozhtsov, Numerical Solutions for Partial Dif- 
ferential Equations, CRC Press, New York, 1996. 

V. G. Ganzha and E. V. Vorozhtsov, Computer-aided analysis of difference 
schemes for partial differential equations, John Wiley & Sons, New York, 
1996. 

Gil' M. I., Stability of finite and infinite dimensional systems, Kluwer Aca- 
demic Publishers, Boston, 1998. 

Gil' M. I., Difference Equations in Normed Spaces, Stability and Oscilla- 
tions, Elsevier, Amsterdam, 2007. 

Edwige Godlewski and Pierre- Arnaud Raviart, Numerical Approximation 
of Hyperbolic Systems of Conservation Laws, Springer- Verlag, New York, 
1996. 



45 



[22] S. K. Godunov, A finite difference method for the numerical computation 
and discontinuous solutions of fluid dynamics, Mat. Sb., 47 (1959), pp. 
271-306 (in Russian). 

[23] L. Gosse, A Well-Balanced Flux- Vector Splitting Scheme Designed for Hy- 
perbolic Systems of Conservation Laws with Source Terms, Computers and 
Mathematics with Apphcations 39 (2000) 135-159 

[24] B. Hanouzet and R. Natalini, Global Existence of Smooth Solutions for 
Partially Dissipative Hyperbolic Systems with a Convex Entropy, Arch. 
Rational Mech. Anal. 169 (2003) 89-117. 

[25] A. Harten, J. M. Hyman, and P. D. Lax;, with appendix by B. Keyfiz, On 
finite-difference approximations and entropy conditions for shocks. Comm. 
Pure Appl. Math., 29 (1976), pp. 297-322. 

[26] A. Harten, High resolution schemes for hyperbolic conservation laws, J. 
Comput. Phys., V. 49 (1983), pp. 357-393. 

[27] Ami Harten, On a class of high resolution total-variaton-stable finite dif- 
ference schemes, SIAM J. Numer. Anal., Vol. 21, No. 1 (1984), pp. 1-23. 

[28] Ami Harten, Uniformly High Order Accurate Essentially Non-oscillatory 
Schemes, 111, J. Comput. Phys., V. 71 (1987), pp. 231-303. 

[29] Robert Horvath, On the monotonicity conservation in numerical solutions 
of the heat equation, Appl. Numer. Math., 42 (2002), pp. 189-199. 

[30] Shi Jin, Runge-Kutta Methods for Hyperbolic Conservation Laws with Stiff 
Relaxation Terms, J. Comp. Phys. 122 (1995), 51-67. 

[31] Shi Jin, Lorenzo Pareschi, and Gioseppe Toscani, Uniformly accurate dif- 
fusive schemes for multiscale transport equations, SIAM J. Numer. Anal. 
Vol. 38, No. 3 (2000), pp. 913-936. 

[32] Shi Jin and C. David Levermore, Numerical Schemes for Hyperbolic Con- 
servation Laws with Stiff Relaxation Terms, J. Comp. Phys. 126, 449-467 
(1996). 

[33] David Kahaner, Cleve Moler, and Stephen Nash, Numerical methods and 
software, Prentice-Hall, New Jersey, 1989. 

[34] N. N. Kalitkin, Numerical Methods, Nauka, Moscow, 1978 (in Russian). 

[35] L.M. Kocic and G.V. Milovanovic, Shape Preserving Approximations by 
Polynomials and Splines, Computers Math. Applic. Vol. 33, No. 11, pp. 
59-97, 1997. 

[36] A. G. Kulikovskh , N. V. Pogorelov, A. Yu. Semenov, Mathematical prob- 
lems of numerical solution of hyperbolic systems, Nauka, Moscow, 2001 (in 
Russian). 



46 



[37] Alexander Kurganov and Eitan Tadmor, New High- Resolution Central 
Schemes for Nonlinear Conservation Laws and Convection-Diffusion Equa- 
tions, Journal of Computational Physics 160, 241-282 (2000) 

[38] Alexander Kurganov and Doron Levy, A third-order semidiscrete central 
scheme for conservation laws and convection-diffusion equations, SIAM J. 
Sci. Comput., Vol. 22, No. 4 (2000), pp. 1461-1488 

[39] P. Lancaster, Theory of Matrices, Academic Press, New York, 1969. 

[40] Randall J. LcVcquc, Finite volum,e m,ethods for hyperbolic problems, Cam- 
bridge University Press, Cambridge, 2002. 

[41] Xu-Dong Liu and Eitan Tadmor, Third order nonoscillatory central scheme 
for hyperbolic conservation laws, Numer. Math. (1998) 79: 397-425. 

[42] L. A. Monthe, A study of splitting scheme for hyperbolic conservation laws. 
Journal of Computational and Applied Mathematics 137 (2001) 1-12. 

[43] L. Mirsky, An Introduction to Linear Algebra, Dover Publications, New 
York, 1990. 

[44] K. W. Morton, Numerical Solution of Convection-Diffusion Problems, 
Chapman & Hall, London, 1996. 

[45] K. W. Morton, Discretization of unsteady hyperbolic conservation laws, 
SIAM J. Numer. Anal., Vol. 39, No. 5 (2001), pp. 1556-1597. 

[46] Giovanni Naldi and Lorenzo Pareschi, Numerical schemes for hyperbolic 
systems of conservation laws with stiff diffusive relaxation, SIAM J. Numer 
Anal., VoL 37, No. 4 (2000), pp. 1246-1270. 

[47] Greg F. Natcrcr and Joso A. Cambcros, Entropy Based Design and Analysis 
of Fluids Engineering Systems, ), Taylor and Francis Group, Boca Raton, 
USA, 2008. 

[48] Haim Nessyahu and Eitan Tadmor, Non-oscillatory Central Differencing 
for Hyperbolic Conservation Laws, Journal of Computational Physics, Vol. 
87, No 2., April 1990, pp. 408-463. 

[49] J. M. Ortega and W. C. Rheinboldt, Iterative solution of nonlinear equa- 
tions in several variables, Academic Press, New York, London, 1970. 

[50] V. V. Ostapenko, On the strong monotonicity of nonlinear difference 
schemes, Comput. Math. Math. Phys., 38 (1998), pp. 1119-1133. 

[51] Lorenzo Pareschi, Central differencing based numerical schemes for hyper- 
bolic conservation laws with relaxation terms, SIAM J. Numer. Anal., Vol. 
39, No. 4 (2001), pp. 1395-1417. 



47 



[52] Lorenzo Parcschi and Giovanni Russo, Implicit-Explicit Rungc-Kutta 
Schemes and Applications to Hyperbolic Systems with Relaxation, Journal 
of Scientific Computing, Vol. 25, Nos. 1/2, November 2005, pp. 129-155. 

[53] Lorenzo Pareschi, Gabriella Puppo, and Giovanni Russo, Central Runge- 
Kutta Schemes for conservation laws, SIAM J. Sci. Comput., Vol. 26, No. 

3 (2005), pp. 979-999. 

[54] Richard B. Pcmbcr, Numerical methods for hyperbolic conservation laws 
with stiff relaxation I. Spurious solutions, SIAM J. Appl. Math., Vol. 53, 
No. 5, pp. 1293-1330, October 1993. 

[55] Richard B. Pember, Numerical methods for hyperbolic conservation laws 
with stiff relaxation II. Higher-order Godunov methods, SIAM J. Sci. Com- 
put., Vol. 14, No. 4, pp. 826-859, July 1993. 

[56] William H. Press, Brian P. Flannery, Saul A. Teukolsky, William T. Vetter- 
ling. Numerical Recipes in C, The Art of Scientific Computing, Cambridge 
University Press, New York, 1988. 

[57] R. D. Richtmyer and K. W. Morton, Difference Methods for Initial- Value 
Problems, 2nd edn, Wiley-Interscience, New York, 1967. 

[58] A. A. Samarskiy, On the monotone difference schemes for elliptical and 
parabolic equations in the case of not self-conjugate elliptical operator, Zh. 
Vychisl. Mat. Mat. Fiz., 5 (1965), pp. 548-551 (in Russian). 

[59] A. A. Samarskii, The theory of difference schemes, Marcel Dekker, New 
York, 2001. 

[60] A. A. Samarskiy and A.V. Gulin, Stability of Finite Difference Schemes, 
Nauka, Moscow, 1973 (in Russian). 

[61] Susana Serna and Antonio Marquina, Capturing shock waves in inelastic 
granular gases. Journal of Computational Physics 209 (2004) 787-795. 

[62] J. H. Wilkinson, The algebraic eigenvalue problem, Oxford University Press, 
London, 1969. 

[63] Huiqing Xie, Hua Dai, On the sensitivity of multiple eigenvalues of non- 
symmetric matrix pencils, Linear Algebra and its Applications 374 (2003) 
143-158. 



48 



